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ABSTRACT 


Providing a simplified model of real terrain has applications to route planning for 
robotic vehicles and militarv maneuvers. In this thesis I explore planar-patch surface 
modeling to represent terrain in a simple and effective way. In planar-patch surface 
modeling the terrain is subdivided into a set of planar subregions. The homogeneity of 
the gradient within a planar subregion simplifies calculating the cost of traversing the 
region, thus simplifving route planning. I have explored three main strategies to model 
the surface: Joint top-down and and bottom-up, strict bottom-up, and presmoothing 
bottom-up approaches. Results of the algorithms are shown graphically by using the 


APL and Grafstat packages, verifying their correctness and accuracy. 
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I. INTRODUCTION 


A. BACKGROUND 

There are many ways to represent a three-dimensional object, usually involving 
complicated mathematical and graphical theories. Representing and displaying a three- 
dimensional object like a terrain is important in applications such as cartography, route 
planning for robotic vehicles, road construction, and planning of military maneuvers. 
A simplified model of real terrain helps the user better understand and analyze it. 

The simplest approximations are piecewise planar, and the very simplest is a 
polvhedron with triangular faces. Assume that a terrain 1s approximated bv à two- 
dimensional array of evenly spaced elevations. Each group of three adjacent points in 
the array defines a planar polyhedral face, and a set of triplets constitutes a planar-patch 
approximation of the surface. Unfortunately, this approach creates many small triangu- 
lar regions, which can make the searching process of route planning too complicated. 

In this thesis, I developed three algorithms for planar-patch modeling of evenly 
spaced elevation data, and each has its own advantages and disadvantages. I have used 
the least-squares method to fit a plane over the points of a terrain, a quadtree algorithm 
to subdivide regions, a gradient clustering method to build regions, and smoothing al- 
gorithms on elevation data to improve continuity of adjoining regions. 

One major advantage of planar-patch modeling over other three-dimensional surface 
modeling is that terrain within a patch is homogeneous in the sense of a constant gra- 
dient. This homogeneity simplifies calculating the cost or speed of traversing the region, 
thus making it easier to analyze for route planning. 

Several thresholds affect my algorithms, and different threshold values will result in 
different models. These include the standard deviation of a fitted plane from given 
points, the ratio of magnitude between adjacent points, and the root mean square dif- 
ference of coefficients of planar equations between adjacent regions. The final result of 
this thesis will be provided to Ron Ross (Ph.D. student, Computer Science Department, 


Naval Postgraduate School) for his research in route planning for robotic vehicles. 


B. ORGANIZATION 
Chapter 2 introduces terrain modeling techniques, including least-squares planar 


fitting and cell classification. 


Chapter 3 is the main chapter of this thesis, defining the planar-patch terrain model. 
Joint top-down and bottom-up and strict bottom-up planar-patch terrain modeling 
Strategies are described in detail. Programming techniques such as the quadtree method 
used in the joint top-down and bottom-up algorithm and the region-growing method 
used in the strict bottom-up algorithm are explained. Finally, the issue of continuity of 
adjacent planar-patches along their common boundaries 1s discussed. Attempts to solve 
this problem, such as smoothing of original elevation data and grouping of the same 
tvpe of points into one region are explained. 

In Chapter 4 comparisons are made between the two algorithms in terms of accu- 
racy and simplicity of terrain representation. Special features and the limitations inher- 
ent to each algorithm are discussed. Experimental results of joint top-down and 
bottom-up, strict bottom-up, and presmoothing bottom-up algorithms are presented in 
graphical form in this chapter. Pictures of the modeled terrain are shown in three di- 
mensions using pictures drawn by the Grafstat and APL packages. Pictures of the 
boundaries of subregions are also included. 

Finally, Chapter 5 summarizes the contributions made by this research and discusses 
some of the possible areas for further research based on this work. The Pascal source 


program is in the Appendix. 


H. REVIEW OF TERRAIN MODELING 


The most common method of terrain representation is by a set of functional ex- 


pressions such as 
fey) = eee 3x y — 5y + y 


where f(x,y) is the elevation function for terrain and 


x and y are the latitude and longitude coordinates. 


This permits the simple calculation of the elevation given any combination of x and y. 
It is known that anv continuous surface can be approximated with an arbitrarily small 
error by a polynomial of sufficiently high degree. Such polynomials can be derived by 
least-squares surface-fitting methods. 

Surfaces mav be defined in tabular form where the value of f 1s given for a selection 
of representative arguments. To find an arbitrary value of f, a table lookup can be per- 
formed to determine the value of the nearest known points and use them to approximate 
the desired value by interpolation. [Ref. 1: p. 212] 

It is possible to specify a surface in terms of guiding points by a generalization of 
the Bezier polynomials or the B-splines. The Bezier polynomial or B-spline method finds 


an approximate curve that passes near a set of given points. [Ref. 2: p. 247] 


A. TERRAIN CELL CLASSIFICATION 

I. QUADRATIC APPROXIMATION OF SURFACE BY LEAST SQUARES 
FIT 
Fitting a least-squares quadratic to a 3 by 3 square of points on a surface is done in the 


following way. The fit as a function of x and v is: [Ref. 3: p. 55] 
f(x) 9 &, 4 ox + kay + kaxó + koxy + key” 


The square of the error between the fit equation and the given terrain points is: 


E’ = X + kox + kay + kax ksxy + key” — flaw)’ 
× y 


The partial derivatives of the squared error with respect to each constant are: 
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To minimize the squared error, each of the above equations is set to zero and simplified. 
Those equations can be further simplified bv assuming that the coordinates of each of 
the points in the 3 by 3 matrix, except the center point, is + one distance unit from the 
center point. By the above assumption, any terms in above equations which contain a 
summation of an x or y coordinate with an odd power will go to zero. Simplified 


- 


equations are: 
0=k DD +k x6 + kç x 6— DD Jx) (11) 
X y Xx y 


0- kx 64» 9 foxy)x (12) 
x y 


Qm dd feo (13) 
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O=ksx4—) > Axyxy (15) 
× y 


0= ۸, × 6 + k4 × 4 + ۸ x6— ) P fag (16) 
x y 
From above equations, if we solve for k,, k,, A;, k,, k; and k,, we get a quadratic 
equation fitted to the 3 by 3 square of points. 
2. TERRAIN CELL CLASSIFICATION BASED ON QUADRATIC 
APPROXIMATIONS 
The eigenvalues of the Hessian matrix for the quadratic function approximated 


fare 3 by 3 points on the surface are: [Ref. 3: p. 61] 


, Ip DK o KOKE) 
m D a e e 


Using the signs of the eigenvalues, we can classify the central point bv Table 1 


on page 6. [Ref. 3: p. 69] 


Table 1. CLASSIFICATION OF SURFACE 
POINTS BY EIGENVALUES OF THE 
HESSIAN MATRIX 


sign of /, sign of 2, 
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B. PLANE APPROXIMATION OF A SURFACE BY LEAST SQUARES FIT 


The elevation z of any point (x,y) can be expressed as: 
z = fix) 
The plane that we want to fit over the surface is: 
2 = ax + by + c 


The square of the error between the fitted plane and the real terrain 1s: 


E= X Y (ax + by + e — fo) 
ay 


Ihe partial derivatives of the squared error with respect to a,b and c are: 





cE = py De +by+c -Ao h> (21) 
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lo minimize the squared error, each of the above equations is expanded and then the 


partials are set to zero. 
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Here, we have three unknowns, a,b and c, and three equations. So we can solve for 


them. 


C. THE GRADIENT AT A POINT 

One of the primary criterion we will use for grouping points into a region is the 
gradient. The gradient of a point is a vector that has a magnitude (slope) and a direction 
for its components. We can represent a three by three elevation matrix of terrain points 
as in Table 2 on page 8, where x and v are the coordinates of a point on a terrain and 


f(x,y) is the elevation of the point. 


Table 2. A MATRIX OF TERRAIN 
POINTS 


f(x-1,v T 1) RSD TS 


f(x t L,y) 


ET 1,3-1) 





The gradient vector of a point (x,y) can be approximated in three ways. The first 


method 1s: 
A, =Ax + 1,5) - Axy) 
A) =fxy + 1)-Axy) 


magnitude(x y) = m + Aj) 
0-6 E 
direction(x,y) = tan CA ) 
1 
The second method is: 
A, — f(x t dy) — fix — ly) 


A, =flxy + 1) y>) 


The magnitude(x,y) and direction(x,y) are the same as above. 


The third way of approximating gradient is by first fitting a quadratic to the three 


by three square of points and then using the coefficients of the quadratic : [Ref. 3: p. 60] 


magnitude(x y) = A a 0 


direction(x y) = tan. ( UE ) 
2 
where k, and k, are the coefficients of the quadratic described in section A. 
The first method considers only one quadrant of the 3 by 3 matrix, thus 1s biased. 
The second method uses a point from each of the quadrants, and is better than the first 
method in the sense that it utilizes more unbiased information. However, it gives an 


mcorrect value in the case as shown in Table 3. 


Table 3. AN EXAMPLE OF A TER- 
RAIN POINTS MATRIX 





When we calculate the gradient of center point using the second method, it will give 
the magnitude value zero and direction value 45 degree, which is wrong. The third 
method does not have such problem, since it utilizes all eight neighboring points. In this 
case it gives the magnitude value 0.2357 and direction value 45 degree. I used the third 


method in my program. 


The gradient vector of a terrain point points in the x-y direction of maximum slope 
of that point; magnitude/x,y) is the maximum slope value and direction/x,y) is the di- 


rection of that slope in the x-y plane. 
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IM. PLANAR-PATCH TERRAIN MODELS 


A planar-patch terrain model is a model that represents terrain in the form of a 
polyhedron. This model is different from the common planar terrain model used in 
graphics in that size of a planar-patch is not fixed, and the shape of a patch does not 


have to be regular. 


A. IMPLEMENTATION 

I chose Pascal as an implementation language for this research because of my fa- 
muliarity and experience, and the Waterloo Pascal on the IBM 370/3033AP mainframe 
was the Pascal compiler that was used. 

The APL language and the Grafstat graphics software package Were used to verify 
the correctness and accuracy of the algorithms. Even though I chose Pascal as my im- 
plementation language, the matrix form of the elevation data suggests APL. With the 
Grafstat software package which runs in the APL environment, one can check the 
intermediate results of an APL program by displaying pictures. Since I experimented 
with different threshold values in my algorithm, APL and Grafstat were useful. The 
three-dimensional pictures in chapter 4 were done with Grafstat. 

Two arrays of records are used as main data structures bv all mv algorithms: 
point array and subregion array. The point, array is a two-dimensional array which 
corresponds to the points in the real terrain. Each element of the arrav 1s a record which 
has the elevation, the subregion ID, and the magnitude of the terrain point as its fields. 
From the values of elevation and magnitude, the subregion ID is determined and as- 
signed to the point. The subregion_array is a one-dimensional arrav of records which 1s 
created during the subregion creation phase. The indices of this array are the ID’s of 
subregions. It keeps necessary information about the subregions such as the equation 
of the plane fitted to the subregion, the boundary of the subregion, and other subsidiarv 
information. 

The general structure of my first two modeling algorithms is shown in Figure 1 on 


page 12. 
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Figure 1. Block Diagram of the First Two Program Approaches 


B. JOINT TOP-DOWN AND BOTTOM-UP PLANAR-PATCH TERRAIN 
NIODELING 

My first algorithm consists of two parts. The first part is the top-down subdivision 
of the original elevation data matrix into subregions and the second part is the 
bottom-up merging of adjoining subregions when adjoining subregions have similar 


plane equations. 


l. TOP-DOWN PHASE 

A quadtree region subdivision method 1s used during the top-down phase of the 
firtst algorithm. The method generates a quadtree by recursively dividing a two- 
dimensional region into quadrants. Each node in the quadtree has four data elements, 
one for each of the quadrants in the region. [Ref. 4: p. 215] The original region is given 
some standard tests such as comparison to a threshold of the maximum magnitude or 
the minimum magnitude of the points within it, depending on the apphcation. Only if 
a region fails the tests and is not of minimum size, it is further subdivided into four su- 
bregions. If a region has both length and width four cells or more, it is not of the min- 
imum size. See Figure 4 on page 16 for the different kinds of minimal subregions. The 
subdivision process described above is repeated for each subregion. Figure 2 on page 
14 shows graphically this algorithm. 

The code works whether a region has a dimensions of odd or even numbers. For 
example, consider 4 by 4 and 5 by 5 matrix regions to subdivide; Figure 4 on page 16 
shows how this method works. 

2. THRESHOLDS USED 

Three thresholds are used in the top-down phase: low_bound, upper_bound and 
standard_deviation (SD). 

The threshold upper_bound is set by the maximum chmbing capability of the 
vehicle in cases where the polvhedral model will be used for route planning. If the ab- 
solute value of the minimum slope of a region is greater than the upper_bound, which 
means every cell has a slope greater than the upper. bound, then the region 1s assumed 
being untraversable bv the vehicle, and the region need not be subdivided further be- 
cause a good polyhedral fit just doesn’t matter for that region. The threshold low. bound 
means that when the absolute value of the maximum slope of a region 1s less than the 
low_bound, which means every cell has a slope less than the low_bound, then the slopes 
in the region are unimportant. Such regions need not be subdivided either. 

For the region which has maximum slope greater than the low. bound and mini- 
mum slope less than the upper. bound, the region is subdivided in four parts according 
to the quadtree algorithm to better 1solate the unlevel slopes. Then we fit a plane to the 
region and calculate the standard deviation of the plane from the real terrain points. 
The standard deviation (SD) is: 


e Y Y LEN oc hus o 
X y 
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QUADTREE REGION SUBDIVISION METHOD 


AN EXAMPLE OF SUBREGIONS CREATED BY THE QUADTREE METHOD 


Figure 2. Quadtree Region Subdivision Method 


PROCEDURE MAKE_SUBREGIONS(YL,YU,XL,XU : INTEGER; 
prev_region : region_ptr); 


(X YL : LOWER BOUND IN Y AXIS, YU : UPPER BOUND IN Y AXIS, 
XL : LOHER BOUND IN X AXIS, XU : UPPER BOUND IN X AXIS X) 


VAR 
F : REGION_PTR; 


BEGIN 
tila > 2) AND ((XU-XL) > 2) THEN CX NOT MINIMUM SIZE X) 


if (pd.max_mag_of_region > low_bound) then (X test X) 
begin 
make_subresions((Cyu-yl) div 2)tyltl , 
CTU xl div 2 )+xlt lee xu  p2; (X Ist quadrant X) 
make subregions(yl , ((yu-yl) div 2)tyl, 
GGxu-xl) div 2otxltl 7 xu , P); (X 2nd quadrant X) 
make subregions(yl , ((yu-yl) div 2)tyl, 
x<1 , اير ع‎ ALO giv 2)txl, p); (X 3rd quadrant X) 
make_subregions(((Cyu-yl) div 2)+y1+1 , yu, 
XI Cu oO diy 2)txl., Pp); (x Gth quadrant X) 


p := prev_region ; (X goto parent region X) 
end 
else stop subdivide (X region passed standard test X) 
else stop subdivide (X region has minimum size X) 


END 





Figure 3. Pascal Code for the Quadtree Region Subdivision Method 


ere /(.x,1') is the elevation of point (x,¥) 
ax + by +c 15S-the plane equation of the region 


nis the total number of points in that region. 


The region is subdivided only if the calculated SD is bigger than the threshold 
SD and it does not have the minimum size. 
3. BOTTOM-UP PHASE 
We now merge adjoining regions having similar plane expressions. Ihe reason 
is that in the process of applying the quadtree subdivision, we often subdivide a region 
which should not have been so completely subdivided. For example, sav a region has 
very steep slopes in the first quadrant and the rest of the quadrants are flat. Since the 


original region will fail in the maximum slope test due to the slopes in the first quadrant, 





a en 


DIFFERENT KINDS OF MINIMUM SIZE 
SUBREGIONS 
| 





BOUNDARIES 





Figure d. Miinimum-Size Subregions Created by the Quadtree Method: This pic- 


ture shows quadtree subdivision method applied to the regions which 
have dimensions of odd or even numbers. and also shows boundaries 


and different kinds of minimum size subregions. 


the top-down phase will subdivide it in four. However, the first, second and third 
quadrants are all flat and should not be separate. 

The adjoining regions of a region are found from the point, array bv looking up 
the subregion ID field of adjoining point. The bottom-up phase 1s done in two sub- 
phases: row_combine and col_combine. The row_combine finds and merges adjoining re- 
gions Which are located in a same row, and the col_comnine finds and merges the regions 


in a Same column. 


Figure 5. Subregions after the Nlerging Process: These subregions are created 
bv performing the merging process over the initial subregions created bv 


the quadtree algorithm. 


To merge adjoining subregions which have similar plane equations, a value 
equation_difference 1s calculated from the two plane equations and compared to the 


threshold abc_difference. The equation_difference ts: 


y 


equation_difference = la — a) 4 (b — by + (Oc) 


P 


where a, b, and c, are the coefficients of plane, and 
az. b, and c; are the coefficients of plane,. 
If this 1s less than the threshold then the two adjoining regions pass the test. 


Then we fit a plane over the two subregions and calculate the standard deviation of the 
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plane, and only if it is less than the threshold standard_deviation we do merge the two 


subregions. 


C. STRICT BOTIOM-UP PLANAR-PATCH TERRAIN MODELING 
l. POINT-CLASSIFICATION PHASE 

In strict bottom-up modeling, we first classify every point in a region bv at- 
taching two tags to it. The first tag classifies slope of the point: level, safe-slope or 
unsafe-slope. A threshold low bound defines the low limit of slope and another threshold 
upper. bound defines the upper limit. The second tag classifies the curvature of the point, 
by the terrain cell classification criteria described in chapter 2. The designators are hill, 
depression, plane and special. The tag special refers to all the rest of the cell classifica- 
tions except the above three. A threshold called curvature is compared to the eigenvalue 
of the Hessian matrix; if its absolute value is less than the threshold, then the eigenvalue 
is regarded as zero.[Ref. 3: p. 69] We group the contuguous points with the same first 
and second tag to form subregions. 

Next the thresholds magnitude ratio and abc, difference are used to combine the 
subregions. —. The magnitude ratio, is used in the region-growing phase and the 
abc difference 1s used as in joint top-down and mottom-up modeling. 

2. REGION-GROWING PHASE 

Ihe region-growing phase consists of two subphases: one-dimensional and 

two-dimensional. 
a. ONE-DIMENSIONAL REGION-GROWING 

The one-dimensional phase starts with the first point (pj) in the first row 
of the point array. The point becomes the head of a linked list and the subregion ID 
field of the point is given the value one. The subregion ID starts from one and is in- 
cremented every time a new subregion is created. Whenever a new subregion is created, 
the necessary information about the subregion, such as the coordinates of the starting 
point and boolean value that savs whether the subregion is active or not, is stored in the 
subregion array using the subregion ID as an index. Then it examines the next point ( 
Di) in the same row. If the ratio computed using magnitudes of gradient ( 
abs(magnitude, — magnitude, ١ 
DEC DU HEUTE REED 1s less than the threshold magnitude ratio, then they 
are linked together to form a linear subregion (R,,). That is, the next, point field of the 
first point in the point_array is given the x and y coordinate of the second point and the 
subregion_ID field of the second point is given the same value as the subregion ID of the 


first point. The average magnitude of the gradient of the points within the subregion is 





Figure 6. One-Dimensional Region-Growing 


Seeulaced and is stored in the ave mag field of the subregion record in the 
Saureson array. [hen repeatedly for each next point (p,,) in the same row, the ratio 
between the average magnitude of the gradient of the subregion and the magnitude of 
the gradient of the new point 1s compared to magnitude ratio. If on this basis new point 
(p,,) should not join the subregion, a new region starts from that new point and the new 
point becomes a head of another linked hist. This process continues for all rows. See 
Figure 6 for an example of this algorithm. 
b. TWO-DIMENSIONAL REGION-GROWING 

The two-dimensional subphase starts with the first linear subregion (R) 
Which was created by the one-dimensional subphase. Actually, it starts from the first 
point in the first column of the point_array. which belongs to the subregion one. It ex- 
amines the point below in the same column, which belongs to another subregion, to 
decide whether the average magnitude of the gradient of the subregions are similar. The 
information, average magnitudes of the gradient of the two subregions, is retrieved from 
the the subregion, array using the subregion ID of each point. If the ratio 1s less than 


magnitude. ratio then they are linked together. The tail of the first subregion 1s linked 
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Figure 7. Two-Dimensional Region-Growine: This algorithm is applied after the 


one-dimensional region-growing algorithm. 


to the head of the second subregion, the tail of the second subregion is left nil, all the 
subregion IDs of the points within the second subregion are changed to the first subre- 
gion ID, and the second subregion in the subregion_array 1s marked inactive. When 
those two linear subregions are merged. the resulting subregion is not one-dimensional 
any more, except if both subregions were single points. Then repeatedly for each next 
point in the same columin, the ratio between the average magnitudes of the gradient of 
the subregions 1s compared to the magnitude_ratio. If a new subregion does not join the 
old subregion. the control is moved to the new subregion and the process starts from 
that new subregion. This process continues for all columns. This is illustrated in 


Figure 7. After the two-dimensional phase, we fit a plane to each subregion using 


least-squares. Then the threshold abc_difference is used to merge adjoining regions that 
have similar plane expressions. The merging is done in two steps. The first step is to 
scan each row of the point, array looking for adjacent points in the same row which have 
different subregion IDs. For every pair Of adjacent subregions, the plane expressions 
are consulted from the subregion_array by using the ID numbers of the subregions as 
the indices. The equation_difference of the two plane equations are calculated and 
compared to abc, difference. If it less, the two subregions are merged. point of the sec- 
ond subregion to the tail point of the first subregion. The second step is identical except 


that scanning for merging is done by column. 


D. AN ALGORITHM IMPROVING THE CONTINUITY OF ADJOINING 
REGIONS 

One common problem with planar-patch terrain modeling 1s that planar patches for 
two adjacent subregions seldom yield a surface that is continuous along the line seg- 
ments bisecting the points used to create the regions. Since the initial requirements for 
this research rule out any other terrain modeling than planar-patch, we had to come up 
with some idea to ensure continuity of the patches. One simple solution 1s to model 
terrain entirely with small triangles, each of which exactly fits three adjacent terrain 
points, as described in Chapter 1. The number of triangles necessary is very large. 

But a minor modification to the two modeling approaches explored in this research 
can give smoothness. One can create supplementarv triangular subregions to fill the gap 
between adjoining subregions. The gap between adjoining subregions have two lines, one 
for each side of the gap. The lines are the edges of the adjacent planes. We can define 
two triangles on the gap bv grouping the first end point of the first edge to the second 


edge and second point of the second edge to the first edge. See Figure 8 on page 22. 
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TWO TRIANGLES CONNECTING 
ADJACENT SUBREGIONS 


SUBREGION 1 SUBREGION 2 





Figure $. Filling the Boundary with Two Triangles 


This is better than the many-triangle model, but is still not optimal, because the number 
of subregions will be increased considerably by the supplementary triangular subregions. 
We have made a different attempt to solve this problem using smoothing of the el- 
evation arrav and an expansion to three dimensions of Rolle's theorem. The basic idea 
ıs that if we fit planar-patches to an evervwhere-convex or evervwhere-concave region, 
the fitted patches are much more likely to be continuous. This means we should trv to 
ensure that all the points within large regions are classified as either planar points, hill 
points or depression points. See Figure 9 on page 23. Smoothing the elevation data 
will help: other kinds of points will tend to be changed into one of the three. 
Unfortunately, because of time, we Were not able to solve ts rather compre, 


problem completely. It is recommended for future research. 


Do 


du A 


INTERSECTION 


THE INTERSECTION TENDS TO LIE BETWEEN THE TWO REGIONS 
IF WE FIT PLANES OVER THE EVERY-WHERE CONVEX OR 
EVERY-WHERE CONCAVE TERRAIN 


Figure 9. How Surface Convexitv Helps 


IV. DISCUSSION 


A. ACCURACY ANALYSIS 

The joint top-down and bottom-up and Strict bottom-up approaches produce fairly 
good terrain models in terms of accuracy and simplicity. Here accuracy means whether 
the model represented the real terrain without losing anything significant, and simplicity 
means the number of subregions. 

The overall accuracv of a model can be verified by comparing the pictures of the real 
terrain and model. But there is another thing that describes the accuracv of a model, 
Standard deviation. We can evaluate the local accuracy in terms of the standard devi- 
ation of the elevation difference for a point in the real terrain and the corresponding 
point in the model. In the top-down approach, each subregion must have a fit standard 
deviation less than the the threshold standard_deviation. Let’s assume that the elevation 
difference between a point on the real terrain and a point on the model has a normal 
probability distribution. If we know the standard deviation of the elevation differences 


of a subregion: 


1. About 68 percent of the points will have elevation difference less than plus or mi- 
nus one standard deviation. 


to 


. About 98 percent of the points will have elevation difference less than plus or nu- 
nus two standard deviations. 


ریا 


Practicallv all points (99.73 percent) 
will have elevation difference less than plus or minus three standard deviations. 


So for example. if we limit the model to have the worst case local elevation differ- 
ence, say 10 feet, with 98 percent certainty, we can set the standard deviation threshold 
to 5 feet. Consequently everv subregion will have a standard deviation less than 5 feet, 
and 98 percent of the points will have an elevation difference less than 10 feet. 

The threshold standard deviation was not used in the strict bottom-up approach so 


the above argument does not hold. 


B. LIMITATIONS OF THE MODEL 

In the joint top-down and bottom-up approach, the minimum number of points in 
a subregion is four and the shape of initial subregion must be a rectangle. This re- 
Striction is due to the quadtree. These limitations prohibit the approach from picking 


up linear (one-cell-wide) terrain features and creating linear subregions. This limitation, 
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however, does not exist in the strict bottom-up approach owing to the linked list of 
points used to represent a region. When a linear region is a straight horizontal or ver- 


tical line, one cannot calculate a plane equation. 


C. EXPERIMENTAL RESULTS 

We now present the output of our terrain modeling program in graphical form. It is 
applied to sample elevation data acquired from the Fort Hunter Liggett terrain data 
base. Shown are the pictures that are produced by our modeling programs using differ- 
ent thresholds. 

As was mentioned in chapter 3, in joint top-down and bottom-up modeling. the 
thresholds were low_bound and upper_bound for gradient magnitude, standard_deviation, 
and abc_difference. In strict bottom-up modeling. the thresholds were low_bound and 
upper bound of gradient magnitude, curvature, abc difference, and magnitude ratio. 

Three sets of pictures are given sequentially in the following section. The first set 
is produced by the joint top-down and bottom-up approach (referred as the first ap- 
proach in the following pictures), the next by the strict bottom-up approach (referred 
as the second approach), and the last by the bottom-up approach applied to the 
smoothed data (referred as the third approach). 

In the array pictures, each subregion has a unique code. If a number does not have 
a prefix then the number ts actually the subregion ID. If it does have a prefix then the 


number has to be converted. See the conversion example below. 


23 ==> subregion ID 23. 
23 = => subregion ID 123. 
23 ==> subregion ID 223. 
-23 ==> subregion ID 323. 
QO ==> either an outermost boundary 


or inside area of a region. 
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EXPERIMENTAL RESULTS 1 OF IHE FIRST APPROACH 
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This is the real terrain that were used in these 


Figure 10. Sample Real Terrain: 


experiments. 
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standard_deviation= 3. The above picture shows the initial subregions 


created by quadtree algorithm. 
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Subregions Created after the Merging in the First Approach: Thresh- 


olds: abc difference 10; standard deviation 2 3; Some of the initial 


subregions created by quadtree algorithm are merged by merging 


process, thus vielding less subregions. 
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Figure 14. 


Terrain Model Created by the First Approach Showing Boundaries 
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Subregions Created by the Quadtree Method in the First 
Approach: Thresholds: low_bound=0.1; upper bound=0.6: 
standard_deviation=7. The above picture shows initial subregions cre- 


ated by the quadtree algorithm using a different threshold (the 
standard_deviation is changed). 
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Figure 18. 


Terrain Model Created by the First Approach Showing 
Boundaries: The planar patches are different from the planar patches 


produced from the previous experiment. 
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Subregions Created by the Region Growing Algorithm in the Second 
Thresholds: 


maenitude_ratio=0.1. This picture is produced after performing the 


Figure 19. 


Approach: low bound-O0.1; upper bound-O.6; 


one-dimensional region growing process (row region growing). 
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Figure 21. Terrain Model Created by the Second Approach: Final product of the 


second approach. 
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Figure 22, 


Terrain Model Created by the Second Approach Showing Boundaries 
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Figure 23. Subregions Created by the Region Growing Algorithm in the Second 
Approach: Thresholds: low bound-0.l; | upper bound- 0.6; 
magnitude ratio —0.2. This picture is produced after performing one- 
dimensional region growing using a different threshold (the 
magnitude ratio is changed). 
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Figure 24. Subregions Created after the Merging in the Second 


Approach: Threshold values used: abc_difference= 20. Two dimen- 
sional region growing is applied and then merging 1s performed using 


different thresholds (abc difference and magnitude_ration are 


changed). 
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Figure 25. Terrain Model Created by the Second Appraach: Final product of the 
second approach. Notice that this picture is different from the picture 


produced by the previous experiment. 
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Figure 26. 


Terrain Model Created by the Second Approach Showing Boundaries 
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5. EXPERIMENTAL RESULTS OF THE THIRD APPROACH 
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Figure 27. Smoothed Sample Terrain: The original sample terrain is smoothed. 
Bottom-up modeling will be performed on this smoothed terrain. 
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Figure 28. Subregions Created by the Third Approach: 


0 0 0 


Thresholds: 


low bound - 0.1: upper bound —- 0.6; magnitude ratio 7 0.1. 
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29. Subregions Created after the Merging in the Third 


Approach: Threshold: abc_difference= 10. 
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Figure 30. Terrain Model Created by the Third Approach: This picture looks 
almost the same as the the picture produced by the second approach; 


this approach produced less subregions than the second approach. 
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Figure 31. Terrain Model Created by the Third Approach Showing Boundaries 
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V. CONCLUSIONS 


In this research, we explored three types of planar-patch terrain modeling. The purpose 
of such models is a simple and yet accurate view of terrain. Since our models were de- 
veloped for route planning, we required a minimal number of subregions, since this di- 
rectly affects the cost of the route planning. 

The first algorithm (joint top-down and bottom-up algorithm) takes the longest ex- 
ecution time among the three. The second (strict bottom-up) and the third (bottom-up 
algorithm applied to the smoothed elevation data) algorithms take almost same time, 
though the third algorithm requires the input elevation data to be smoothed first, re- 
quiring some eXtra preprocessing time. As far as result accuracy (simplicity) is con- 
cerned, the first algorithm 1s better than the others. The second and third do not use a 
standard deviation threshold, resulting in possible inaccuracy, and tend to produce more 
subregions. Since all of the three algorithms use almost the same data structure, a two- 
dimensional point_array and a one-dimensional subregion_array, they take almost the 
same space. As for the application environments of these algorithms, the first algorithm 
could be used in situations where the accuracy is critical, while the second and third 
algoritms could be preferred in situations where a short execution time is required and 
linear terrain features should be identified. 

Since the algorithms used in this research operate on large arravs and involve heavy 
computations, a multiple-processor computer architecture would vield superior process- 
ing capability since the same operations could be performed for each part of the input 
matrix simumtaneously. As I conclude this research, I regret that I could not solve the 


continuity issue completely. 
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APPENDIX A. PASCAL CODE OF JOINT TOP-DOWN AND 


BOTTOM-UP ALGORITHMI 


(*$s350000*) 
( seared sede sede re ede jede de dede eese ye jede e dede vede jede deve TENE dee ye VEN dede vede Te ye ye dede ede jede ee de eee ede ye ye ee eee eee 


Author = Yee, Seung Hee. 
Date = 3 June 1988. 
Input = 1. Elevation data acquired from Fort. Hunter Liggett 
data base (40 by 40 square). 
2. file of magnitude of gradient of the above elevation 
data created by using the program in Appendix 3. 
Output = 1. Elevation data of the model created by this program. 


2. Elevation data of the boundary points, other points 
are given the value one so that when we display 
this model in three-dimensional pictures, we can 
only see the boundaries of the patches. 

Computer = Waterloo Pascal on IBM 370/3033A4P, 
Naval Postgraduate School. 
Description = 
This program is the Pascal implementation of the 


Joint top-down and bottom-up terrain modeling. 
WIM IER TEIN ye yeso ye dede ye je ye ye vede eve ye ede e dee jede ye ve ve ye veo de vede eed eve jede dede eye ede dede dede jede jede jeje dese ©) 


program top, down(input,output); 
const 
murnft 


- 808; 
low bound 


2 
(ue (* mag > 10 % slope then divide *) 
(* mag <= 10 % then disregard *) 
upper bound - 0.6; (* mag < 60 % then do not subdivide *) 
max_sd = 3; 
abc_difference =10 ; 
max_group_num= 300; 


3: 


row = 40; 
col = 40; 
increment = 3; 


type 
elrange = 1.. 10000; 
pointrec = record 
el : elrange ; 
mag: real; 
gp: OT max group num; 
end; 


magrec = record 
mag : real; 
LC lI TOW, 
end; 


point, array-array (.1..row,l..col.) of pointrec; 


grouppointer = -grouprec ; 
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grouprec = record 
LIL EM 
el. cr al... row: 
magmax,magmin : magrec ; 
up = Srouppointen, 
end; 
constarray = array (.1..3,1..4.) of real; 


nodepointer = -node; 
node = record 
r,C : integer; 
next : nodepointer; 
end; 


coefrec = record 
erl cerh, 
ECUSCEL O MASON 
a,b,c : real ; 
minmag,maxmag : real; 
sd : real ; 
ee,ex,ey : real; 
side : 0..max_group_num; 
constarr : constarray; 
active : boolean; 
list : nodepointer ; 

end; 


subregion_array = array (.1..max_group_num. ) of coefrec ; 


var 
points point array: 
groups : grouprec ; 


subregion : subregion_array ; 
ok : constarray; 
data,datam,doutl,dout2 : text; 
counter,r,c,i : integer; 
avg,ee,ex,ey : real; 


top,p : grouppointer; 
temparr: array (.1..max group num.) of boolean; 


function variance (a,b,c,ee,ex,ey: real; ok: constarray): real; 


begin 
variance := (ee-(2*a*ex)-(2*b*ey)-(2*c*ok(.3,4.)) 
+(ara*ok(.1,1.))+(2*a*b*ok(.2,1.)) 
+(2*b*c™ok(. 2,3. J)+(2*a™c™ok(. 1,3. )) 
+(b*b*ok(. 2,2. ))+(c*c*ok(. 3,3. ))) 
¡CORSO 
end; 


procedure initialize; 


var 
r,c :integer; 
begin 
counter := 0 ; 
IOC CSI COCO OO 
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end; 


proce 


var 


601 


for r 
begin 


end; 


dure 


1 to row do 
pointsl.r.)( cJ- £p 


:= 0; 


3 


= 1 to max group, num do 
subregion(.r.).a := -999; 
subrecion(t- r8 59:2 7-999: 
subregion(.r.).c := -999; 
subregion(.r.).active := false; 
subregion(.r. ).side := 0; 
subregion(.r.).sd := -999; 


(* procedure initialize *) 


getconst(var points: point_array; 
var arr: constarray; 


i, j: integer; 


rena eler Integer; 
var ee,ex,ey :real); 


begin 
oT 1 := 1 to 3 do 
begin 
for j) -= 1 -to 4 do 
arrr, j): =Ù; 
end; 
ee := 0; 
ex := 0; 
ey := 0; 
moe 1 := rl to rh do 
begin 
on) = CLIO CI do 
begin 
lo) arr 1,1 +H i*i): 
I DG a RTD 
arn 13. J: QATI. L3: + 1; 
AE Arr 1 4 i oomnts(.r.].).e1*1); 
ATE OL amp 2, iJ 2*1); 
AGM we wen —4Ep092 2.)50]*1): 
AE A ATI 2) j 
are + a (polmts(.1,)].).el*3); 
a irr(.3,1.)9535 
ange 3.2. ): zartr(33,2. )t 1: 
SHE 3,3. Jy sarr(. 9,93. )*1; 
ee -u-adspecc4scpointso 12].).el; 
ee := ee + sar(points(.i,j.).el); 
exe — eX teGpoints(-1,].).el * 1); 
ey := ey + (points(.i,j.).el * j ); 
end (* j *) 
end; (+ 1 7) 
end; (* procedure get const *) 


E 


procedure gauss(ok: constarray; var aa,bb,cc: real); 


var 
1,j : integer; 
temp : array (.1l..4. ) of real; 
t : real; (* error checking purpose *) 
begin 
qoum 
while ok(. 1,1.) 2 0.0 do i := 1i # 1 ; (F&F TOW rotation ~) 
for j:= 1 to 4 do 
begin 
temp(. j.) := ok(.1,3.); 
ACA O AO 
okl i) S tenp. j). 
end; 
t := ok(. 1,4); 
forsi = 2 tomode 
for j := 4 downto 1 do 
ok(.i,j.) += ok(. 1. ]- SokC dI غ4‎ ec eee 
i := 2; 
while ok(.1,2.) = 0.0 do i1i := i + 1; 
for j := 1 to 4 do 
begin 
Lemp...) = OKI 
OKC OR ea ee 
ôk- ici) S tenp ار‎ 
end; 
for j := 4 downto 1 do 


ok(.3,j.) := ok(.3,3.)*ok(.2,2.) + ok(. 2, j. )*(-ok(. 3,2. )); 


if ok(..3,3. )/ <> 0.0 then 


begin 
ok(.3,45.) v2 0kC-29, 0 IS OEIL 
ok(..2,4.) := ("ok(.2,45.) S OK on USER MO O 
ok(.1,4.) := ( ok(.1,4.) - ok(.2,4. )*ok(. 1,2.) 
- ok(s99 56 )sokRC ISI TORO S 
aa := ok(.1,4.); (* solution **) 
bb := 0k(.2,4.); (FT solution *) 
ce := ok(-3,4.);. (= solucions) 
(* gauss error = t-( aa*ok(.1,1. )+bb*ok(.1,2. )+ cc*ok(.1,3.)) *) 
end; 
end; (* procedure gauss *) 


procedure prefix(num: integer); 
begin 
case Fy div 100) mod 10) of 

write(' 2 

Wriel TF 
write(':'); 

write( -'): 
write('+'); 


Ono 
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5 rite D 
6 : write('*'); 
7 write( & ): 
Saa i 
9 : write('z'); 
end; 
end; (* procedure prefix *) 


procedure printsubl( points: point_array); 
var 
EF C,i,j,k : integer; 
begin 
page; 
for k := 1 to counter do 
with subregion(.k.) do 
tore c—erltl to crh* bdo 
tor j -= CcltI to cer- l do 
Points- i,j: ): gp- = 0; 
forr := 1 to row do 
Bor ec := 1 to col do 
begin 
if c<> col then begin 
BECTIX(PpOINES(. TC. J SP); 
were (Points r e SBP MOG 1000:2); 
end 
else begin 
BECI POINTS r,e) 8P); 
wriceln (points. r,c.)-gp mod 100): 2); 
writeln; 
end; 
end; 
end; (* procedure printsubl *) 


procedure printsub2(points:point, array); 
var 

ENS corl;r2,cl,c2 : integer; 

current : nodepointer; 
begin 

for r:= 2 to row-1 do 

for c:= 2 to col-1 do 

poudts( r,c.).gp :- 0 ; 
for k:= 1 to counter do 
if subregion(.k.).active then 


begin 
current :- subregion(.k.).list ; 
while current <> nil do begin 
rl: =current>. x; 
cl: =current>. c; 
if current>.next <> nil then begin 
2 —CUrrene-, NeXt. T; 
c2: =current>. next>. c; 


end 
else begin 
TOM EL: 
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C2, 1-061. 
end; 
if ri 2 r2 then begini: 


r2; r2:-60m gez v. end: 


if cl > c2 then begin ¢:= ¢2:%e2: ="el ie end 


if rl=r2 then for c: cl tome zmac 
points(.rl,c.).gp := k; 

if cl=c2 then for r := 
points( spo ci) sp sso 


TI co rdo 


current := current. next; 
end; 
end; 
for r := 1 to row do 
fore = lI CO CON dao 
begin 


if c<> col then begin 

prefix(points( Ec.) 2p). 
write((points(.r,c.) -gp moa 100) 1 

end 

else begin 
prefix( points( ccm. 
writeln((points(.r,c.).gp mod 100): 2); 
writeln; 

end; 

end; 
end; 
(* procedure printsub2 *) 


procedure get nodes(var subregion: subregion array); 


type 

direction = ( up,down,right,left ); 
var 

1 : integer; 


current,head,p : nodepointer; 
temp : direction ; 


procedure get next(var p : nodepointer; 
var temp : direction); 
var 
i,j : integer; 
done : boolean ; 


procedure do_dirívar done: boolean ; dir:direction ); 


begin 

new(p); | 

poor -É ue 

5[ >: ©6.دم 

done := true; 

case dir of 
up : temp := up; 
down : temp := down; 
right : temp := right; 
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left : temp := left ; 
end; (* case *) 
end; (* sub_sub_ procedure do_dir *) 


ein sub procedure get next *) 


done := false ; 
O. E 
دار‎ BPC; 


case temp of 
right : begin 
igo 
while not done do begin 
pones tia ep” k 
then do_dir(done,up); 
if (not done and 
Gpolmmies (aaa cp) <> ik) ) 
then do dir(done,down); 
i Subregion. k.) lista. r = i) 
and (subregion(.k.). list>.c = j)) then 


begin 
new(p); 
DIETS), 
CE 5 
done := true; 
end; 


if not done then j:= j+1: 
end; (* while *) 
end; 
left : begin 
dle 
while not done do begin 
1ل‎ pOimis@ itl. |. ). ep =k 
then do_dir(done,down); 
il not done and 
üpodntso ci jeg op sk) ) 
then do_dir(done,up); 
if not done then j:= j-1; 
end; (* while *) 
end; 
up : begin 
iS ie 
while not done do begin 
if points(- i ]} l). 2p =k 
then do_dir(done,left); 
if (not done and 
(Points. i=l, j 999p s k)) 
then dol dir done, right); 
if not done then i:= i-1; 
end; (* while *) 
end; 
down : begin 
i := itl; 
while not done do begin 
TERPO n Se iTi) gp Fk 
then do_dir(done,right); 


EID 


if (not done and 
(points(. itl, ji. <> MED) 
then do dir(done,left); 
if not done then i:= itl; 
end; (* while *) 
end; 
end; (* case *) 
end; (* sub procedure get_next *) 


begin (* procedure get nodes”) 
for k := 1 to counter do 
if subregion(.k.). active then 
begin 
new(p); 
p^.r :-subregion(.k.).crl ; 
p*.c :-subb5egion kh) ccu 
while ( points(.p^.t-1l,p^.c.).gp = k) do 
PD E ات له‎ (* move to the upper most row of the gp *) 
while ( points(.p^.r,p^.c-1.). gp = k) do 
Doc Sp cei (* move to the left most col of the gp *) 


subregion(.k. ). list :=p ; 
current := p ; 
temp. = rights; 
repeat 

get_next(p,temp); 

current-^.next :-7 p ; 

current := : 
until (subregion(.k. ). list=.r = p^.r) 

and (subregion(.k. ). list^.c = p¬.c) ; 


p™. next := nil; 
end; 
(* for k := 1 to counter do 
if subregion(.k. ). active then 
begin 
write('gp (| ,k:3,'**==>'); 
current := subregion). k DF listi. 


while current <> nil do begin 
write(current>. 1 2.00 current CIS) 
current := currentn^. next; 
end; 
writeln; 
end”) 
end; (* procedure get_nodes *) 


procedure makesub(rl,rh,cl,cr : integer; prev: grouppointer); 
type 
why = (nondivisible,magnitude,std,steep); 


var 


dif, rC- TInNteger:; 
flag : why; 
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sdoa b.cc real: 


procedure stopdivide(flag : why ); 


var 
r C : integer; 
begin 
counter := counter + 1; 
with subregion(. counter. ) do 
begin 
CTI FEL 
crh := rhe 
cel :s cl; 
Cer = cr. 
maxmag := p”.magmax. mag ; 
minmag := p”.magmin. mag ; 
end; 
for e: iel 0066 
forr =r] to rh do 
pounts(.r.9(. Cc. ). gp >= counter; 
(* case flag of 
nondivisible : 
write('* non divisible area!'); 
magnitude 
write('* stop divide! by mag. J); 
Sta ; 
oee 5 52652 it fits well: 
Secor : 
0215) s stop it is too steep. ); 
end; 


rien: area Il rh cl cr Il. rbi 3/6153. 6r: 35 
has gpnum=' counter: 3 ); 


*) 
p := prev ; 
end; (* sub procedure stop divide *) 


begin (* procedure makesub *) 
new(p); 
p”. magmax. mag := -999; 
p^.magmin.mag :- 999; 
pr] --rl; 
Drm :— rh, 
posue cl; 
BICI = Cr; 
p^. up :7 prev ; 
ror eC. CI lO Cr do 
Or TTI TO fh do 
begin 
if points(.r.)(.c.).mag > p>.magmax. mag then 
begin 
p7.magmax. mag := points(.r. )(.c. ). mag; 
p”.magmax.r := r; 
p”.magmax.c := cj 
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end; 
if points(.r.)(.c.):-mag < p. magmin mae Then 


begin 
p”. magmin. mag := points(.r.)(.c.).mag; 
p^.magmrn. t = r; 
p7. magmin. c := €; 
end; 
end; 


(* max,min magnitude for region rl,rh,cl,cr before division *) 


getconst(points,rl,rh,cl,cr,ok,ee,ex,ey); 
gauss(ok,a,b,cc); 
sd := sqrt(variance(a,b,cc,ee,ex,ey,ok)); 


(* sd for gp before division *) 


if ((rh - rl) »-2) and (Cer - cl) > ZF Then 
if (p^.magmax.mag ? low bound) and 
(p^.magmin.mag < upper bound) then 
if sd > max sd then 
begin 


makesub( rl. , )) 25221( div 2 FETS 
cl , OCcreel) divo2) tcl re 
makesub(((rh-rl) div 2) + rl + 1 , rh, 
Cl, CTCL) OI 2 eer ا‎ 
makesub(rl , ((rh-rl) div 2) + rl, 


((cr-cl)vdrzv 2) cl BUNTE 
makesub(((rh=-nl) diy 2) + 21 olen 
(Ccrscl) div 2) + cl +.) EE ا‎ 
p := prev ; 
end 
else stopdivide(std) 
else stopdivide( magnitude) 
else stopdivide(nondivisible); (* smallest area *) 
end; (* procedure makesub *) 


procedure do_combine(var pt: point_array; 
var subregion: subregion_array; 
var last,new : integer); 


Var iis). 2156 2 
tall,current — 0). max _2roupenum, 
begin 
subregion(. last. ). active := true; 
subregion(.new. ). active := false; 
current := subregion(. last. ).side ; 
if current = O then 
begin 
subregion(. last. ).side := new; 
current := new; 
end 


else begin 
while current <> 0 do 
begin 
tail := current ; 
current := subregion(.current. ).side ; 
(* traverse *) 


58 


end; 


(* subregion arr *) 


subregion(. tail. ).side := new; 
(* to link new gp *) 
current := new; (* into existing link *) 
end; 
while current <> 0 do (* reassigning gp # to new region *) 
begin 
for ii := subregion(. current. ).crl to 
subregion(. current. ).crh do 
Or jj -~= sübregionl. current.) cci to 
subregion(. current. ).ccr do 
PEUS- gp = last; 
current := subregion(. current. ). side; 
end; (* while *) 
Or 1 :-— 0955.3 do 
Hor jj :—- 1 to 4 do 


subregion(. last. ). constarr(.ii,jj.) := 
subregion last. )-constarrl.ii,jj:-) +t 


Sllbrecion(. new. ). constarr(. 11, )}. ); 
subregion(. last. ).ee := coefs(. last. ).ee + 
subregion(. new. ).ee ; 
subregion(. last. ).ex := coefs(. last. ).ex + 
subregion(.new. ). ex ; 
subregion(. last. ).ey := coefs(. last. ).ey + 
subregion(.new. ). ey ; 


with subregion(. last.) do 
begin 
gauss(constarr,a,b,c); 
sd sqrt(variance(a,b,c,ee,ex,ey,constarr)); 
end; 
end; 


(* sub procedure do combine *) 


procedure combine(var pt: point, array; 
var subregion: subregion, array); 
Var 
tabc difference,i,j,last,new : 
joinabc, joinsd, joined : 


integer; 
boolean; 


procedure compareabc(var join: boolean); 


var 
temp : real ; 
begin 
temp := sqrt( sqr(subregion(. new. ).a - 
subregion(. last. ).a) + 
sqr(subregion(.new.).b - 
subregion(. last. ).b) + 
sqr(subregion(.new. ).c - 
subregion(. last. ).c) ); 


if temp <= tabc_difference then join : 


else join := false; 


end; 


crue 


(* sub procedure compareabc *) 


procedure checksd(first,second: coefrec; var join: boolean); 


var 
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ii,jj : integer; 


begin 
for 11, >=. 1 6 Oo 
for jj := 1 to 4 do 
first. constar SÓ 
firsSt.constarr(o 141998 F 
second. constarrí 175 1]5)5 
first.ee := first.ee + second.ee ; 


first.ex := first.ex + second.ex ; 
first.ey := first.ey + second. ey ; 
with first do 
begin 
gauss(constarr,a,b,c); 
sd >= sqrt(variance(a,b,c,ee,ex,ey,constarr)); 
if sd > max_sd then join := false else join := true ; 
end; 
end; (* sub procedure check_sd_before_combine *) 


begin (* procedure combine *) 
tabc_difference := increment; 
while tabc_difference <= abc_difference do begin 
joined := false; 
for i:= 2 to (row-2 edo 
begin 
last := pt 12 eo, 
fOr}: = 2 EO (COI CNIS 
begin 
new :— pto TD PD: 
if ((i2row-2) and (j2col-2)) and (not joined) then 
subregion(.new. ). active := true; 
if (new<>last) and 
(subregion(. new. ). maxmag <= upper_bound) and 
(subregion(. last. ).maxmag <= upper_bound) then 
begin 
compareabc( joinabc); 
1f joinabe 
then begin 
checksd( subregion(.]last.), 
subregion(.new.),joinsd ); 
if joinsd then begin 
do, combine(pt, 
subregion,last,new); 


joined := true; 
end 
else begin 
subregion(. last. ). active := true; 
joined := false ; 
last := new ; 
end; 
end 
else begin 
subregion(. last. ). active := true; 
joined := false ; 
last := new ; 
end; 


end 
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else begin 


Subregion(, vase. ). active += true; 
joined := false ; 
last := new ; 
end; 
end; 
end; (* row combine finished *) 
Bor |] := 2 to (row-2) do 


begin 
last. := pt. 2, 07-9925 5 
Gcr 15 2 to (col 2) do 
begin 
new -= ptl-ifl,j-): 2p; 
if (new<>last) and 
(subregion(.new. ).minmag <= upper_bound) and 
(subregion(. last. ). minmag <= upper_bound) then 
begin 
compareabc( joinabc); 
if 01285 
then begin 
checksd( subregion(. last. ), 
subregion(.new. ), joinsd ); 
if joinsd then begin 
do_combine(pt, 
subregion, last,new); 


joined := true; 
end 
else begin 
subregion(. last. ). active := true; 
joined := false ; 
last := new ; 
end; 
end 
else begin 
subregion(. last. ). active := true; 
joined := false ; 
last := new ; 
end; 
end 
else begin 
subregion(. last. ). active := true; 
joined := false ; 
last := new ; 
end; 
end; 
end; (* column combine finished *) 
tabc_difference := tabc_difference + increment; 


end; (* while *) 
end; (* procedure combine *) 


procedure cal avg dev(var subregion:subregion array; 


var points:point array; var avg:real); 
var 


sum: real; 
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1.1, 2p: integer, 


begin 
sum := 0.0; 
tor 1:= 2 to cow-lede 
for j:= 2 to col ١ 
begin 
gp: poiuts( 15]. o 
sum := sum + abs(( subregion(. gp. ).a*i + 
subregion(. gp. ).b*j + 
subregion(. gp. ).c) 
= points(. 1i, 0 E 
end; 
avg := sum / ((row-2)*(col-2)); 
writeln('average deviation = '*,avg :15); 


end; (* sub procedure cal_average_deviation *) 


procedure writedata( points: point_array); 
var 

ELO: integer: 

tempa,tempb,tempc : real; 


begin 
for r:= 2 to row-1 do 
for C:= 2 to COLL do 
begin 
tempa := subregion(. points(.r,c. ). gp. ). a; 
tempb := subregion(.points(.r,c. ). gp. ). b; 
tempc := subregion(. points(.r,C. ).gP. ).c; 
points(.r,c.).el := trunco (tempat 1 tempb-c tempe): 
end; 
tor c := 1 tO Col Oo 
for r := row downto 1 do 
writeln doutl points r ceo: 
end; 


(* procedure writedata *) 


procedure writetempdata(points:point array); 
Var 
k,r,c,rl.r2,cl.c2 TECO 
tempa,tempb,tempc : real; 
current : nodepointer; 
begin 
IOI T= 2 COTO IAG 
for C= 2 CO Col LEG 
points(.r,c.).el := 1 ; 
for k:= 1 to counter do 
if subregion(.k. ). active then 


begin 
current := subregion(.k. ). list ; 
while current <> nil do begin 
Dl: =current ta. r: 
cl: =current>. c; 
if current>.next <> nil then begin 
r2. current» 2268. و1‎ 


co Current NEXt rac, 


end 

else begin 
r2 se rl: 
CES Cl; 
end; 


then begin TC: 2 T2. r2:= rl; rl:= r; end;‏ 2خ - 211 دل 


AA 6202 Cenn begin C=C? CIS Cl, cl:= c; end; 
subregion(.k. ). a; 
tempb := subregion(.k. ).b; 
tempc := subregion(.k.).c; 
ieri r2 then for c := cl to c2 do 
points lcs -cruncsCctempa^*rl - tempb*c + tempc); 
ii Cle) Chen IOI 1t :— rl to r2 do 
points(.r,cl.).el := trunc (tempa*r + tempb*cl + tempc); 
current.: = Current. next; 
end; 
end; 
for € i= 1 to col do 
Hom © ;= row downto 1 do 
"tryzteln(dout2,points(.r,c.).el); 
end; 
(* procedure writetempdata *) 


cempa : 


procedure check; 
var 
mum, r,C,i : integer; 
begin 
for c:= 1 to counter do temparr(.c.) := false; (* check *) 
num := 0; 
for i:= 1 to counter do 
if subregion(.i.). active then 
begin 
T= Î; 
temparr(.r.) := true ; 
num := num +1); 
while r <> 0 do 
begin 
with subregion(.r.) do 
temparr(.r.) := true; 
r:= subregion(.r. ).side ; 
end; 
end; 


for c:= 1 to counter do temparr(.c.) := false; (* check *) 
num := 0; 
for i:= 1 to counter do 
if subregion(.i.). active then 
begin 
w= i; 
temparr(.r.) := true ; 
num := num +1 ; 
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while r <> 0 do 
begin 
with subregion(.r.) do 
temparr(.r.) := true; 
r:= subregion(.r. ). side ; 
end; 
end; 


writeln('num of gps = ',counter :3,' * lowbound=' ,low_bound: 3, 
' upper_bound=" ,upper_bound:3,' * max_sd=',max_sd:3 ); 
writeln('unexamined gps are ==>'); 


c:= 0; 

for i:= 1 to counter do if temparr(.i.) = true then c:= ctl 
else write('*', 1:3. *'): 

writeln; 


writeln('# of gps examined =>',c:4,' * abc_difference=', 
abc_difference); 
writeln('# of combined gps =>',num :4); 


end; (* procedure check *) 


( N === === ===“ == === === === == =>========= 


gag OO ا‎ 2222225 22222 2 a aa aea $ $ ) 
begin 
page; 
Maatialize ; 
new(top); 
Hops. up := nil ; 


Meset(data. data ell): 

reset(datam, data magl'); 

rewrite(doutl,'data outl'); (* data after combine process *) 
rewrite(dout2,'tdata outl'); (* data of boundary *) 


for c := 1 to col do 

for © := row downto 1 do 
readin(data,points(.r.)(.c.).e1); 

Lor c:- 1 to col do 

for r := row downto 1 do 


rFeadin(daram.points(.r.)(-c.).mag); 


mákesub(2,(row-1),2,(col-1),top); 
printsubl(points); (* print the subregions created by makesub *) 


for i := 1 to counter do 
with subregion(.i.) do 
begin 


getconst(points,crl,crh,ccl,ccr,ok,ee,ex,ey); 
gauss(ok,a,b,c); 


sd = sqrt(variance(a,b,c,ee,ex,ey,ok)); 
(* sd for region i after division *) 
constarr := ok ; 
end; 


combine( points ,subregion); 
get_nodes(subregion); 
printsub2(points); (* print the subregions created by combine *) 
writedata(points) ; 

(* write modified elevation data into file ‘data out#'*) 
writetempdata(points); 

(* write elevation data into file ‘tempdata GUÊ ~) 
cal avg dev(subregion,points,avg); 
(* check; *) 


end. 
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APPENDIX B. PASCAL CODE OF STRICT BOTTOM-UP ALGORITHM 
(*$s450000*) 


( Vevezesevevevevevedevevededeseve fede edes eVededevedeyede deze eee vede dededese defe dede fededeveyedeseve dese deve sese dedeyevededeyevedededevo 


Author = Yee, Seung Hee. 
Date = 3 June 1988. 
Input = 1. Elevation data acquired from Fort. Hunter Liggett 
data base (40 by 40 square). 
2. file of magnitude of gradient of the above elevation 
data created by using the program in Appendix 3. 
3. file of curvature of points created by using the 
program in Appendix 3. 
Output = 1. Elevation data of the model created by this program. 
2. Elevation data of the boundary points, other points 
are given the value one so that when we display 
this model in three-dimensional pictures, we can 
only see the boundaries of the patches. 
Computer = Waterloo Pascal on IBM 370/3033AP, 


Naval Postgraduate School. 

Description = 
This program is the Pascal implementation of the 
Joint top-down and bottom-up terrain modeling. 


PRO e‏ اک هاه هاه هاه هله مله له مله داه دک ملو ها "fa elo et an a‏ ماه م اه ھا هله د ل داه د له د له د له هھ اه داه د له د لے م له ذه م أ م كه له م له م لهم لوه أده کے وات وأ ها tratar rl‏ هه اد هاه دص 
A TOv)‏ رم ل م أنه هله هل *4 er ae ae r A ES > E Û ss > EE > e» rr ev‏ مك مم و we LA ʻe er‏ ذم ذه > > je és Y cf: ee + e nt au ee see's eu ed E rS we > ya ri e * y 5: y s Jv es $ vs‏ 


PROGRAM Bottom, up(input,output); 
const 


low bound =5 0.1; (* magnitude < 10 % then level slope “*) 
(* magnitude > 50 % then unsafe slope *) 

up. bound = 0.6; ( else safe slope 30 

magnitude_ratio = 0.2; (* ratio of magnitude : 


dif of two mags / old mag *) 
abc_comparison_const =20; 
max_group_num = 400; 


row = 40; 
col = 40; 
type 
elrange = 1..10000; 
mag type 7 (level,safe,unsafe); 
point = record 
r,c: integer; 
end; 
nodepointer = -node; 
node = record 
r,c : integer; 
next : nodepointer; 
end; 
point_rec = record 
el : elrange ; 
mág : real; 
magtype: mag_type; 
cuc: char: 
£p : O..max group num; 
penext :wpormu 
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end; 

|poumtearray - array (.1:..x0w,l..col.) of point, rec; 
polynomial array = array (.1..3,1..4.) of real; 
coef_record = record 

start : point; 

upper_left : point; 

a,b,c : real ; 

sd : real ; 

ee,ex,ey : real; 

side : 0..max_group_num; 

poly array : polynomial array; 

active : boolean; 

list nodepolnter ” 

numpts : integer; 

avgmag: real; 

end; 
coef array = array (.1..max_group_num. ) of coef_record ; 
var 

points : point_array ; 
subregion : coef_array ; 
ok : polynomial_array; 
data,datam,datac,doutl,dout2,dout3,dout4 : text; 
ccunter,r,c,1,num : integer; 
tempmag : real; 


function variance (a,b,c,ee,ex,ey: real; ok: polynomial_array): real; 


begin 
variance :7 (ee-(2*a*ex)-(2*b*ey)-(2*c*ok(.3,4.)) 
*(a*a*ok(.1,1.))*(2*a*b*0ok(.2,1. )) 
L2 bsctokio 252 ))502*a*c*ok(C. 1,3 )) 
tebsb OKC. 2,2. ) )+(c*c*ok(. 3,3. ))) 
¡CORSO 
end; 
procedure initialize; 
var 
mc : integer; 
begin 
counter := 0°: 
tore => 1: to col do 
forr := 1 to row do 
begin 
Poimes(s rye: ). 2 = 0; 
podnest.r,c. ). AM e 
Points- rc. ).ptnext:.c := 0 ; 
end; 
for r := 1 to max_group_num do 
begin 
subregion(.r.).a := -999; 
subregion(.r.).b := -999; 
subregion(.r.).c := -999; 
subregion(.r.).active := true; 
subregion(.r. ). side := 0; 
subregion(.r. ).sd := -999; 
end; 
end; (* procedure initialize *) 
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procedure getconst( r,c: integer; 
var arr:polynomial array; var ee,ex,ey : real); 
begin 
arr(.1,1.)9?sSargr( L1. Ge). 
arr(. 1,2: ): Sarr- 1,2 tes. 
áàrr(.1,9,)arr( 1.9 ee 
arr(.1,4,): =arr(. 1,4. )+(polmts (ecc 
arr(.2,1.):*—arr(.2 I P ODE 
arr(.2, 2.) =arr( 2 200) 
arr(.2,3.)"=2ar rt. 2 OS 
arr(.2,4.):-arz(.2,5. pointes Er. como 
Arr( go] STE TT 
arro 3 2 —qrp dm DE 
arr(.3,39.): -8rr( 43.5 91m 
arr(.3,4. ):=a2rr(.3,6: )^porsnts( orc eee 
ee :— ee + sqr(points(.r,c. ).el); 
ex := ex t (points(.r,c.)-el 4r); 
ey := ey t (posnesQ rco) vel "cm 
end; (* procedure get const *) 


procedure gauss(ok: polynomial_array; var aa,bb,cc: real); 


var 
153]; “integer: 
temp : array (.1..4.) of real; 
t : real; (* error checking purpose *) 
begin 
1:= 1]; 
while ok(.i,1.) = 0.0 do i:= i3+ 1; (* row rotation *) 
for j} := 1 to 4 do 
begin 
templ. J-) 3 = Ont. ine 
ole Lycos qu 
ORCI 1.) E CENDI 
end; 
€: 3 ="0k(..1,4.); 
for 1s = 2 to 3 do 
for j := 4 downto 1 do 
ok(. i,j.) := ok(.1,3.)*ok(C. 1; lage OR oO apr 
iste 2: 
while ok(.i1,2.) » 0.0 do i:-i t1; 
for j := 1 to 4 do 
begin 
templ. J-) -Sokl 2 sir 
ok(.2,j.) SOR T; 
ok(. i,j.) := temp(. j. ); 
end; 
for j := 4 downto 1 do 
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ok(.3,j-) := 0k(.3,3.)*ok(.2,2.) * ok(.2,3. )*( -ok(. 3,2. )); 


if ok(.3,3.) <> 0.0 then 


begin 
eno) COR) OK. y 
A O 4.) = OK. 3,4. )*ok(.2,3.)) / okt. 2,2. ); 
Oba c ok(.1,4.) - ok(.2,4. )*ok(. 1,2. ) 
oko E32 )) “OR. Ll.) 5 
aa := ok(.1,4.); (* solution *) 
bb += OK. 2,4. ); (* solution *) 
CC = OKC. 3,4.); (Er solution *) 


(*writeln('gauss error =', 
E Aa ol bb ok. 1,2. )F ccřok(. 1,3.)))*) 
end; 
end; (* procedure gauss *) 


procedure prefix(num: integer); 
begin 

case E div 100) mod 10) of 
Write a2 
0 | D. 
wg ted B 
write(!-' 
write('+' T 
write n, 
write( * ); 
write é€ J; 
write. PE 
Write, = je 


0O 0 Os UI EF ني مر نم ين‎ 


end; 
end; (* procedure prefix *) 


procedure printsub(points: point_array); 


var 
Ec integer; 
begin 
page; 
Tol 1 := l to row do 
Hre = to col ao 
begin 


if c<> col then begin 
21501-01250225 
"005-0561250 0 2, اتج :(ت‎ 560 100:25 

end 

else begin 
preto ntsc) 2p); 
Wirateln( (points (.1,¢.). eo mod 100): 2); 
writeln; 

end; 

end; 
end; (* procedure printsub *) 


procedure printsub2( points: point_array); 
var 
ez el,cz : integer; 
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current : nodepointer; 
begin 
page, 
for r:= 2 to row- isco 
for c:= 2 to col-1 do 
pointS(. TC oo ١ 
for k:= 1 to counter do 
if subregion(.k.).active then 


begin 
current := subregion(.k. ). list ; 
while current <> nil do 
begin 


rl! =Currentes: =, 
cl:=curtent- Te: 
if current-.next <> nil then begin 
r2:-currentoonextosr: 
62: €-currentsmextc c: 
end 
else begin 
r2 dde 
cz Mes 
end; 
if ri > T2 then 5682111 1 


= r; end; 
if cl > c2 then begin c:= c2; c2: 


T2: TZ EI NER 
= c; end; 


clic 


Ho 


if rl=r2 then TOI CSCI TO CS 
points (e cio pS 
if cl=c2 then for r := 


points- r, cl ep. kK; 


TI TCO TOG 


current := current”. next; 
end; 
end; 
tor rz lI to row do 
for cee = l to colar co 
begin 


if c<> col then begin 
prefix (points rre.) - 2p); 
write( (points: rc.) -gp mod r100) 2i 

end 

else begin 
prefix points r.c- 2p). 
writeln((points(. r c). gP moa TOO aE 
writeln: 

end; 

end; 
end; 
(* procedure printsub2 *) 


procedure printmagcur; 


var 
r,c : integer; 
begin 
(¥ page”) 
for TL: =|, tO TOW mado 
for cs 1 to coló de 
begin 


70 


if c<> col then begin 
case .ع2 )0125م‎ ). magtype Of 


level urite( 1 poin sAr e). cur; ): 
Sale Write sS 0) Te). cur, ~); 
unsafe :write('u Bpoumcs e c. cur, S 
end; (* case *) 
end 


else begin 
case points(.r,c.).magtype of 
Merem uriteln 1 Smoimts(sc.c.). cur, s 


Sale WrIiltelnl DoSUNPSQ E. C.) GUI,  ); 
unsal ee vrien u pots re. cur, ): 
end; 
writeln; 
end; 
end; 
end; (* procedure printmagcur *) 


mote makeone_rowl(var pointl,point2 : point); 
egin 
Poumes(  pOintl.1,pOlmel.~c.).ptnext.r ¿= point2.r; 
pones pointi. r pomt e). ptnext. c := polnt2.c; 
PODES (PONELE porte e). 2p := 
DOINESI DOIN |i POintiwwe.-). 2p; 
with subregion(. counter. ) do 
numpts := numpts +1; 


end; (* sub procedure make one, l *) 


mere makeone_row2(var pointl,point2 : point); 
egin 
pones (pom pontiac ptuext- r^ - point2.r; 
Poms DOlntl ET polntiaec. ).ptnext.c := pount2.c; 
POEs. DOiNnec. ©,pOInt2. ic. ), coe: = 
Points pote lees pOinel.c. je cm: 
with subregion(. counter. ) do 


begin 
numpts := numpts +1; 
avgmag :— avgmag*(numpts-1)/numpts + 
BONS in point24c.)^maz/numpts; 
end; 
end; (* sub procedure make one_2 *) 


procedure rowlink; 


var 
m : integer; 
5111م‎ 2051222 1 5012 

begin (* sub procedure rowlink *) 
counter := 0; 
for i:= 2 to (row-1) do 
begin 
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counter :2 counter + 1 
points(.1,2.) gp counter: 
subregion(. counter) start 
subregion(. counter. ).start.c : 
subregion(. counter. ). upper_left : 


J 
Dr 


subregion(. counter. ). start ; 


subregion(. 
subregion(. 
Or ji2 
begin 
points, 
point 
Pointe. 
pointa: 


counter. ). avgmag 
counter. ).numpts 
to (col-2) do 


po 
l; 


KOR 
pr je: 


c := j+l; 


ints(. 1,2 mas: 


if ((points(.i,j+1.).magtype <> safe) and 


(points(. i,j. ). magtype 
makeone, rowl(pointl,point2) 


points(. i, j+1. ).magtype) ) then 


(* (level level) or (unsafe unsafe) *) 
safe) and 


else if ((points(.i,j.).magtype 


(points(.i,j+l. ).magtype 


safe)) and 


(points(.1i,J]-)ocur - points dl cuRthen 
if (abs(points(- 1, j+). mace 

subregion(. counter. ). avgmag)/ 

subregion(. counter. ). avgmag) < 


magnitude_ratio then 


makeone, row2(pointl,point2) 


else begin 
counter := counter 
Points- i Eep 
subregion(. counter. 
subregion(.counter. 
subregion(. counter. 
points 
subregion(.counter. 
pointl := point2; 


end 


else begin 
Counter ICONO 
POINCS 1TL EOE 
subregion(.counter 
subregion(.counter. 
subregion(.counter 


+1 ; 


+1 ; 


3 


:= counter; 


DE 
Je 
JE 
e 
Jn 


. ). avgmag : 


start := point2; 
upper_left := point2 ; 
avgmag : > 

Ts mas: 

numpts := 1; 


3 

counter; 
.). start 
). upper, left 


¿= point2; 
:= point2 ; 


points. 1111.2 


subregion(.counter. ). numpts : 


pointl ;:= poinezZ: 
end; 
end; 
end; 


end; (* sub procedure row link *) 


procedure makeone_col(last,new : point); 


var 


iiz JJ- integer; 


12 


1 ; 


paulcurrent : point; 


begin 
subregion(. points(. new. r,new.c. ). gp. ). active := false ; 
with subregion(.points(. last.1r,last.c.). gp. ) do 
begin 


numpts:=numpts + 
subregion(. points(. new. r,new.c. ). gp. ). numpts; 

avgmag: -avgmag*( numpts - 
subregion(.points(.new. r,new. C. ). gp. ). numpts)/ 
numpts + 
subregion(.points(.new. r,new. Cc. ). gp. ). avgmag* 
subregion(. points(.new. r,new.c. ). gp. ). numpts/numpts; 

end; (* reassigning number of points and average magnitude  *) 


if (subregion(. points(.new. r,new.c.).gp.).upper_left.r + 
subregion(. points(.new.r,new.c. ). gp. ).upper_left.c) < 
(subregion(.points(. last. r,last.c. ). gp. ). upper_left.r + 

subregion(.points(. last. r,last.c. ). gp. ). upper_left.c) then 
subregion(.points(. last. r,last.c. ). gp. ). upper_left := 
subregion(.points(. new. r,new. Cc. ). gp. ). upper_left ; 


current := points(. last.r,last.c. ). ptnext ; 
if current.r = 0 then (* or current.c = 0 *) 
begin 


points(.last.r,last.c.).ptnext := 
subregion(. points(.new.r,new.c. ). gp. ).start ; 
current := subregion(. points(.new. r,new.c. ). gp. ).start ; 
end 
else begin 
while current.r <> 0 do 


begin 
tail := current; 
current: == polnts(. current. r,current.c. ). ptnext; 
end; (* traverse ptnext of lastpt link 
to the end of link *) 
pones aca toaiinc. je ptnext:—  (* now link the tail 


points to the head of*) 
subregion(. points(. new. r,new.c. ).gp.).start ; 
(* newpt link *) 
current := (* into existing link *) 
süubregion points new: r new-c. ).gp-).start ; 
end; 
while current.r <> 0 do (* reassigning gp # to new pts *) 
begin 


pones current. current.c.).2p 3= points(. last. r,last.c.). gp; 
current := points(.current.r,current.c. ). ptnext; 
end; (* while *) 
end; (* procedure makeone_col *) 


procedure collink; 


var 
i,j : integer; 
pointi point? : point; 
begin (* sub procedure col link *) 
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flor 1:2 2 tous col. do 


begin 
for 1:7 2 tod rows2)gdo 
begin 
pointi r =mi 
porntl c2: 5 1 
point2. aii 
põint2- cr: =F, 


if (points(.itl,j.).gp <> points Jo (5 رشت‎ 
if ((points(. i+1,j. ).magtype <> safe) and 
(points(.i,j.).magtype = points(. i+1,j. ).magtype)) then 
makeone, col(pointl,point2) 
(* (level,level) or (unsafe,unsafe) *) 
else if ((points(. i,j. ). magtype = safe) and 
(points(. i+1,j.).magtype = safe)) then 
if (points(.i,]. ). cur = points(. it eel 
if abs(subregion(. points(.it+l,j. ). gp. ). avgmag- 
subregion(.points(.i,j.). gp. ). avgmag)/ 
subregion(.points(.i,j.). gp. ). avgmag < 
magnitude_ratio 
then 
maàkeone col(pointl,point2); 
pointl :- point2; 
end; 
end; 


end; (* sub procedure collink *) 


procedure do combine(last,new : point); 


begin 


LOs integer, 

tail, current point, 

subregion(. points(. new. r,new. C. ). gp. ). active := false ; 
with subregion(.points(. last: r. last eker ade 

begin 


numpts: =numpts + 
subregion(. points(. new. r,new.C. ). gp. ). numpts; 
end; (* reassigning number of points *) 


if (subregion(. points(. new. r,new.c. ). gp. ).upper_left.r + 
subregion(. points(.new. r,new.c. ). gp. ).upper_left.c) < 
(subregion(. points(. last.r,last.c. ). gp. ).upper_left.r + 

subregion(.points(. last. r,last.c.). gp. ). upper_left.c) then 
subregion(.points(.last.r,last.c.).gp.). upper left :- 
subregion(.points(.new.r,new.c.).gp. ). upper. left ; 


current := points(. last.r,last.c. ).ptnext ; 
if current.r - 0 then (* or current.c 0 *) 
begin 


points(.last.r,last.c.).ptnext := 
subregion(. points(.new. r,new.c. ). gp. ). start ; 
current := subregion(. points(.new.r,new.c.).gp.).start ; 
end 
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else begin 
while current 1T € UES 


begin 

Tall — current; 

cuüurrent := points current. r, current. c. ). ptnext; 
end; (* traverse ptnext to link the head *) 
points(.tail r, tail c.: ).ptnext:= (* of the *) 


subregion(. points(. new. r,new.c. ). gp. ). start ; 
(* new point *) 
current := (* into existing link *) 
subregion(.points(.new.r,new.c.).gp.).start ; 
end; 
while current.r <> 0 do (* reassigning gp # to new pts *) 
begin 
Ponts current. r, current. c. ).2p : > polnts(. last.r,last.c. ). gp; 
current := points(. current. r,current.c. ). ptnext; 
end; (* while *) 


Bomeii :—- 1 to 3 do 
for jj:= 1 to 4 do 
eueregion(.points(. last.r,last.c. ). gp. ). poly_array(.ii,jj.) := 
Supregion points(. last.r,last.c.). gp: ). poly array(. 11,3]. J+ 
subresion(.polnts(. new. £,new. cc. ). ep. ): poly_array(. ii, jj. ); 
subregion(. points(. last. r, last. c. ). gp. ).ee := 
subregion(.points(. last. r,last.c. ). gp. ).ee + 
subregion(.points(. new. r,new.C. ). gp. ).ee ; 
stibresion(. points(. last. r,last.c. ). gp. ).ex := 
SupEectonml  poines(. last. Y,last.c. ). gp. ).ex + 
subregion(. points(.new. r,new.c. ). gp. ).ex ; 
subregion(. points(. last. r,last.c.).gp.).ey := 
subregion(.points(.last.r,last.c. ). gp. ).ey + 
subregion(. points(.new. r,new.c. ). gp. ).ey ; 
With subregion(. points(.last.r,last.c.).gp.) do 


begin 
sd := sqrt(variance(a,b,c,ee,ex,ey,poly_array) ); 
end; 
end; (* procedure do_combine *) 


procedure combine(var pt: point_array; 
var subregion: coef_array); 
var 
i,j : integer; 
last,new : point; 
joinabc, joinsd, joined : boolean; 


procedure compareabc(var join: boolean); 


var 
temp : real ; 
begin 
temp := sqrt( sqr(subregion(. points(.new. r,new.c. ). gp. ).a 


51121681602). 50125). 1852. 251858. Cc. ). gp. ). a) 
sqr(subregion(.points(.new.r,new.c.). gp. ). b 
subregion(.points(.last.r,last.c.). gp. ). b) 
sqr(subregion(.points(.new.r,new.c.). gp. ).c 


4e‏ و ل ,م 


subregionl. points(. last. Tr laste oo: 
1f temp <= abc_comparison_const then join := true 
else join := false; 


end; (* sub procedure compareabc *) 


begin (* procedure combine *) 
joined := false; 
for 1 := 2 tor row-s2) do 
begin 
Las te 
laste : 
for 4 := 2 to (col?) do 
begin 
new. r := 1; 
new.c := jtl; 
if ((i=row-2) and (j=col-2)) and (not joined) then 
subregion(. points(. new. r,new.c. ). gp. ). active := 
true; 
if (points(.new.r,new.c. ). gp <> points. last r, last cc. joe 
then begin 
compareabc( joinabc); 
if joinabc 


1; 


nou 
N hH. 


then begin 
do, combine(last,new); 
joined := true; 
end 


else begin 
subregion(. points(. last. F,laSt. €. J. 2P. J. active 


true; 
joined := false ; 
last := new ; 
end; (* joinabc *) 


end 
else begin 
subregion(. points(. last. lastic epi) actives 


true: 
joined := false ; 
last := new ; 
end; (* last <> new *) 
end;(* end j *) 
end; (* end i *)(* row combine finished *) 
for j := 2 to (ron =Z) Ao 
begin 
Last. = 2; 
las. - j 
for i:= 2 to (col-2) do 
begin 
new.r := itl; 
new.c := j ; 
if (points(.new:r new.c.). gp ~- points alast r last eeen 
then begin 


compareabc( joinabc); 
if joinabc 
then begin 
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do combine(last,new); 


joined 


end 


else begin 


sübregionl- points. last. r,last. c. ). gp- ).-active : 


:= true; 


joined := false ; 
last := new ; 
end; (* join abc *) 


end 
else begin 


true; 


subregion(.points(. last. r,last.c. ). gp. ). active := 


joined 


last 
end; 


end; 


end; 


end; 


procedure get_boundaries(var subregion: 


typ 


var 


e 


direction = ( up,down,right,left ); 


IR 


current,head,p : 
temp : 


procedure get next(var p: 
Var TEMP. 


var 


(* column combine finished *) 


:= false ; 


:— new ; 


(* last <> new *) 


(* procedure combine *) 


FYE : integer; 


direction ; 


25: integer; 
done : boolean 


procedure do dir(var 


begin 

new(p); 

paro; 

prc r= j; 

done := true; 

case dir of 
up Impe: 
down temp : 
fight e tenp: 
left temp : 


3 


end; (* case *) 
end; (* sub sub procedure do dir *) 


(* sub procedure get next *) 


begin 
done := false ; 
1 


H Hu mW i 


nodepointer; 


coef array); 


nodepointer; 


up; 
down; 


TOP 


left-: 


> 


done: boolean ; 
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direction); 


dir direction ); 


true; 


Jf = Dae, 
case temp of 
right : begin 
و ل‎ 
while not done do begin 
if points(-i-1 M £p 5 
then do_dir(done,up) (* go up and done *) 
else if( (subregion(.k.).list-.r = i) 
and (subregion(.k.).list-.c = j)) then 


begin (* return to Org point 
new(p); 
PT TE 
مر‎ NS 
done := true; 
end 
else if (points(.i,jtl.).gp = k) 
then j:= 3+1 (x og right *) 


else if (points(.itl,}j. ). gp=k ) 
then do_dir(done,down) (* go down and done *) 
else do_dir(done, left); (* go back and done *) 
end; (* while *) 
end; 
left : begin 
while not done do begin 
if points dis) = K 
then do dir(done,down) (* go down and done*) 
else if (points(.1,)- eco — لعا‎ 
then j := j-1 (* go left *) 
else if (points(27-1.)\2)- cp. — 5 
then do_dirídone,up) (* go up and done *) 
else do_dir(done,right); (* go back and done *) 
end; (* while *) 
end; 
up : begin 
15 
while not done do begin 
11 5601256: ترع:(1,351‎ > 1 
then do_dir(done,left) (* go left and done *) 
else if (points(.i-1,j.).gp = k) 
then i := i-l (* go up *) 
else if (points(.i,jtl.).gp =k ) 
then do_dir(done,right)(* go right and done*) 
else do_dir(done,down); (* go down and done *) 
end; (* while *) 
end; 


while not done do begin 
if points (© RK 
then do_dir(done,right)(* go right and done*) 
else if points- itl I) IT kK 
then i := itl (* go down *) 
else if points(.i,j-l.).gp =k 
then do dir(done,left)(* go left and done*) 
else do dir(done,up); (* go up and done *) 
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end; (* while *) 

end; 

end; (* case *) 
end; (* sub procedure get next *) 


begin (* procedure get boundaries *) 
for k := 1 to counter do 
if subregion(.k. ). active then 
begin 
new(p); 
pur supreedgon(.Kk.).upper-left.r ; 
p^.c :-subregion(.k.).upper left.c ; 
hile ( points.. p~ r-l, pi c. ).gp = kK) do 
PE: BINTI, (* move to the upper most row of the gp *) 
02112 ) 256031215) Pp. reps. c-i- ). 9p = k) do 
أن يدم‎ supe ce]; (* move to the left most col of the gp *) 


subregion(.k. ). list := p ; 
011 20511115). peyr,p 7, cr. ). gp = k then 
begin (* starting toward right groups *) 
cürrent -= p; 
temp := right ; 
repeat 
get_next(p,temp); 
current”. next := p ; 
current := ; 
uncii Esubrecion(.K. ). 11st.r = pa.r) 
and (subregion(.k.)>).list-.c = p”.c); 
if (temp=down) and (points(.p”.r+1,p”.c.).gp=k) then 


begin (* right and down points gp *) 
current :=p ; (* hinged on start point *) 
temp := down ; 
repeat 
get_next(p,temp); 
current^.next :-7 p ; 
cürrenl DE 


üncil (subregion. K. ). list r = p=.r) 
and (subregion k: ). lista e = p.c)? 


pi next := nil; 
end 
else 
BE next = hil, CF right but not down gp *) 
end 
else if points(.p”.r+1,p”.c.).gp = k then 
begin (* starting toward down group *) 
current :=p ; 
temp := down ; 
repeat 
get_next(p, temp); 
current-.next :=p ; 
current SDS 


until SubresiOn(. kay list=. r = pas) 
and subregion TK o lista = p.c); 
p^.next :7 nil; 


19 


end 
else 
p*.next :- nii. (*» onec posnuntuepe) 
end; 


( 7 o o ماي‎ o eC EE WE GE GEMENS MEE GERE EE ME 








the following segment of code may be used for checking purpose 
for k := 1 to counter do 
if subregion(.k. ). active then 
begin 
write('Boundary of group ',k:3, is ==>’); 
current := subregion(.k. ). list; 
while current <> nil do 
begin 
write(current 4, ri2, , Weurrentoc. 2. == )- 
current := current>. next; 
end; 
writeln; 
end; 











end; (* procedure get_boundaries *) 


procedure calc_coeffs; 
var 
k,i,]).f.cgpnum proum . simcece:. 
ty tl pos: 
rowchanged,colchanged: boolean; 
begin 
gpnum := 0; 
ptnum := 0; 
for k := 1 to counter do 
if subregion(.k.).active then 
begin 
rowchanged := false; 
colchanged := false; 
fOr Ss 1 fto h? do 
for j:= 1 to 4 do 
subregion(.k. >). poly_array(.i,j.) := 0; 
subregion(.k. ).ee : 
subregión k. ex... 


3 
0; 
subregion(. k. ey m 
ptnum := ptnum +subregion(.k. ).numpts ; 
(* num of linked pt *) 

gpnum := gpnum +1 ; (* num of active gp *) 
t := subregion(.k. ). start; 
While t-r <> O0 36 
begin 

TI STE 

with subregion(.k.) do 

getconst( t.r, t.c,poly array,ee,ex,eyj; 
t: 7 points( -t r.t- C) PENCET. 
if (t.rs>0) and (Eesti en 


rowchanged := true; 
if (t.¢<>0) and (£. cel cc) Sennen 
colchanged := true; 


end; 
with subregion(.k.) do 
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begin 
if rowchanged and colchanged then 


gauss(poly_array,a,b,c) (* normal points *) 
else if not(rowchanged) and not(colchanged) then 
begin 

a := 0; 

b:=0; 

c := poly_array(. 3,4. ); 
end (* one point gp *) 
else if not(rowchanged) and colchanged then 
Dore (* horizontal line group *) 

an 


i j:= 4 downto 2 do 
o ON Je 
poly_array(. 3 »j- Ji poly arraáay(.2,2.) + 
pol mama OT PO Ly array(:3,2.)); 
poly_array(.3,4.) := 
poly_array(. 3 .4.( / 
2501 0 


poly_array(.2,4.) := (poly_array(. 2,4.) - 
poly_array(. 3,4.) * 
poly array 2D / 
poly array 2,2. ); 
b := poly_array(.2,4.); 
c := poly_array(. 3,4. ); 
end 
else if rowchanged and not(colchanged) then 
begin (* vertical line group *) 
b:= 0; 


for j:= 4 downto 1 do 
poly_array(.3,].) := 


poly- array 3j- )7 
poly_array(.1,1.) + 
poly_array(.1,j. )* 
(pola array(. .3,1.)); 


poly _array(- 3,4.) : 
pn 5 / 
BOL ATTA), 
poly_array(.1,4.) := (poly_array(. 1,4. )- 
poly_array(.3,4. )* 
Poly array- 1,3.))/ 
poly array- 1,1. ); 


a := poly array(.1,4.); 
c := poly_array(.3,4. ); 
end; 
(* oa LO B10, e LO) *) 
end; (* with subregion(.k. ) 7e) 
end; (* if subregion(.k.).active **) 


writeln('active pts=',ptnum:5,' * active gpnum =' ,gpnum: 3); 
end; (* procedure calc_coeffs *) 


procedure write planar data(points: point, array); 
var 
TIC’ integer; 
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tempa,tempb,tempc : real; 


begin 

for r:= 2 to row-1 do 

for c:=2 to col= ts 

begin 
tempa := subregion(. points(.r,c. ). gp. ). a; 
tempb := subregion points o r e maA. 
tempe := subregion points -r,e EDE 
points(.r,c.).el := trunc (tempa*r + tempb*c + tempc); 

end; 

for c := 1 to col do 

for r := row downto 1 do 


writeln(doutl poes roc Mou 
end; 
(* procedure write_planar_data *) 


procedure write, boundary. data(points: point, array); 
Var 
k,r,c,rl,r2;cloc2  Wa3ntecem. 
tempa,tempb,tempc : real; 
current : nodepointer; 
begin 
for r:= 2 to row-1 do 
for c:= 2 tO col] lease 
points(.r,c.).el := 1; 
for k:= 1 to counter do 
if subregion(.k. ). active then 


begin 
current := subregion(.k. ). list ; 
while current <> nil do 
begin 


rl: =current”.r; 

cl: =current^. c; 

if current”-.next <> nil then begin 
r2: =current >. next>. 4x3 
c2: =current>. next>. c; 


end 
else begin 
r2 :- rl; 
c2 := cl; 
end; 
if rl > r2 then begin r:= 12722) — ieee = ena 
if cl > c2 then begin c:= c2; c2:= cl; cl:= cj end; 
tempa := subregion(.k.).a; 
tempb := subregion(.k. ).b; 
tempc := subregion(.k. ).c; 
if rl=r2 then for c := cl to c2 do 
points(.rl,c.).el := trunc (tempa*rl + tempb*c + tempc); 
if cl2c2 then for r := rl to 22006 
points(.r,cl.).el := trunc (tempa*r + tempb*cl + tempc); 
current := current”. next; 
end; 
end; 
for C SLES COL do 
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tor rm: = tow downto 1 do 
ecen doutz points. r.c) el); 
end; (* write_boundary_data *) 
(* procedure write_boundary_data *) 
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begin 


page; 
initialize ; 


reseut data, data ell 

reset(datam, data magl'); 

reset(datac, data curl’); 

rewrite(dout1,'data out1'); (* data after combine process *) 
rewrite(dout2,'tdata outl'); (* data of boundary *) 


for. c += 1 CO Col do 

for r := row downto 1 do 
readin(data,points(.r,c.).el); 

for e TITS Col o 

for r := row downto 1 do 

begin 


readln(datam,tempmag); 
if tempmag < low_bound then 
points(.r,c. ). magtype := level 
else if tempmag > low_bound then 
points(.r,c. ). magtype := safe 
else points(.r,c. ).magtype := unsafe; 
points(.r,c. ).mag := tempmag; 
end; 
tor. eo =. 1 to comido 
for r := row downto 1 do 
readln(datac,points(.r,c. ). cur); 


(*printmagcur”) (* print magnitude and curvature *) 
row link; (* link points in row and make subregions *) 
printsub(points);(* print subregions created by rowlink *) 
comer (* link points in row and column and make subregions *) 
calc_coeffs; (* calculate plane expressions for the 
subregions created *) 
combine(points ,subregion); 
(* combine adjacent regions if they 
have similar plane expressions *) 
calc_coeffs; (* calculate plane expressions for the 
final subregions *) 
get_boundaries(subregion); 

(* get boundaries of subregions in linked list *) 
printsub2( points); (* print subregions created by combine *) 
write_planar_data(points); (* write out the elevation of the model *) 
write boundary data(points);(* write out the elevation of 

the boundaries of the subregions *) 
end. 
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APPENDIX C. PASCAL CODE FOR SMOOTHING ELEVATION DATA 
AND CALCULATING GRADIENT 


(*Ss300000*%) 


program calc_slope_curv_smooth( input output); 


( veveseveseveveve deveve eveveve vevevevevevede vede jeveveve Jevedeveve Ve veye ve Vevedevede Je jede deveveveve و‎ dede veveveveve 


Author = Yee, Seung Hee. 
Date = 3 June 1988. 
Input = 1. Elevation data acquired from Fort. Hunter Liggett 
data base (40 by 40 square). 
Output = 1. File of magnitude of gradient of the points. 
2. FIle of curvature of the points (optional). 
3. FIle of smoothed elevation data (optional). 
Computer = Waterloo Pascal on IBM 370/3033AP, 


Naval Postgraduate School. 
Description = 
This program contains two major procedures: 
calc_slope_curvature and smoothing. 
The procedure calc_slope_curvature calculates the slope 
and curvature of each points, the procedure smoothing processes the 
elevation data to get smoothed elevation data. 
Original elevation data is represented in feet and there are 40*40 
data cells in one window. The distance between each cell is12.5 
meters , so it needs to be converted to get a metric system. 


VeTeveveveveveveveve eve vevevevevevevede ve veveveveveve ve ve see jede vedeveveveveveveveve ve ve YedeveYevevedev eve deveveveveeveveve vede deveveve ) 
const 
8 . — m p a 
distance_unit = 1; (03 29089512 5 = 41.010577) 
row = 3; 
col = 3; 


curv_thresh = 0.05; (* threshold for eigenvalues of biquadratic *) 
upbound = 0.6; (* slope > 60 % then unsafe slope *) 


lobound = 0.1; (* slope < 10 % then level slope *) 
type 
elrange = 0..10000; 
pointrec = record 
el : elrange ; 
slope: real; 
cur. char; 
end; 
point_array=array (.1..row,1..col. ) of pointrec ; 
var 
peduts.spoints : point.array ; 
data,outslope,outcur,smoutel,smoutslope,smoutcur > text; 
counter,r,c,i,num: integer; 
k2,k3,k4,k5,k6,slope, laml,lam2: real; 


procedure init(var pt: point_array); 
var 
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Ec, Jj: integer 


begin 
for r := 1 TO TOW OG 
begin 
pt(-r,l.).slopew S5 7 
pt(.r,col. ).stepe., = 99 
pt(. 7,1. )cur 0> 
ptr col cu UE 
end; 
fOr CS I EO Col Oo 
begin 
ptt. 1,6.) slope 1> IO 
DEL. TOW TICA) So pe : > 999 
PEI CIC wea 
pt( stow, caecum == ul 
end; 
end; 


procedure calc_slope_curvature(var pt: point_array); 


var 
r;c,1,]:; integer, 
templ,temp2,termx,termy,dtermx,dtermy, 
templam ,termxx,termyy,termxy,term : real; 
begin 
for r:= 2 to row-1 do 
= 2 


Porc to col-1 do 
begin 
term := 0.0; 
termx := 0.0; 
termy := 0.0; 
termxx := 0.0; 
termyy := 0.0; 
termxy := 0.0; 
for i :=-1 to 1 do 
for J) :=-1 to Tdo 
begin 
term := term + pt(.r+i,c+j.).el; 
termx *= termx + 4peC ricer yee 
termy := termy + ipt TIC ee 


termxx := termxx + j" jJ pt(. rti, ct ed 
termyy := termyy + ixi pt rti a E 


termxy := termxy t ij piel rti CT el. 
end; 
k2 := termx/(6*distance_unit); 
k3 := termy/(6*distance_unit); 
k4 := (termxx/2 - term/3)/distance_unit; 
k5 := termxy/(4*distance_unit) ; 
k6 := (termyy/2 - term/3)/distance_unit; 


laml := ((2*k4+2*k6) + 
sqrt(abs(sqr(-2*k4-2*k6)-4*(2*k4*2*k6-sqr(k5)))))/2; 
lam2 := ((2*k4+2*k6) - 
sqrt(abs(sqr(-2*k4-2*k6)-4*(2*kA4*2*k6-sqr(k5)))2))/2; 
if abs(laml) < abs(lam2) then 
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begin 


templam := laml; 
laml := lam2; 
lam2 := templam; 
end; (* eigenvalues of bigger slopenitude becomes laml *) 


Slopes SGIE(SGEK-2) T Sgr(k3)); 


pt(.r,c.).slope := slope; 
if (laml>curv_thresh) and (lam2>curv_thresh) then 
pate. cur de (* depression  *) 
else if (lami<- cubo -thresh) and (lam2<-curv_thresh) then 
Pre) cur: = h (* hill or peak *) 
else if qu thresh) and 
(abs(lam2)<=curv_ thresh) then 


Pr cO cure = p (* planar *) 
else if e curv_thresh) and (lam2>curv_thresh) then 
PE). cur > 's' (* saddle %) 


else if (laml<-curv_thresh) and 
(abs( lam2)<= =curv _thresh) then 
PRE). cur := r (+ ridge 30 
else if (abs(lam2)<=curv_thresh) and 
(laml>curv_thresh) then 


BECI CUT vy (* valley *) 
else if (lam2<-curv_thresh) and (laml>curv_thresh) then 
CO re e CUr : =. 4 (* pass *) 


end; 


end; 


procedure smoothing(var points,spoints: point, array); 
var 

penis] temp, alí : integer; 
begin 

"oF e S 2 to (row-1) do 


momee 2 to (col-1) do 


begin 
A points. r,e. ). Slope > lobound) 
then begin 
temp := 0; 
TOI 1:= 2-1 to rt! do 
temp :- points(.T3,c. elo temp; 
EHE -c-hlotosckl- do 
temp := points(.r,j.).el + temp; 
spoints(.r,c.).el := trunc(temp/6); 
end 
6156 526313215). 1ه6.زءعءوء‎ := 56125: 2-615 
end; 
for r := 1 to row do 
begin 
speimecuas, |. ).el := points(.r,1.).el; 
Spornts( 1, col. ).el := points(.r,col.).el; 
end; 
Or C= l to col. do 
begin 


SPOILS (le. ).el := points(.1,c.).el; 


87 


= points(.row,c. ).el; 


SpOints( snow cm el 
end; 
end; (* procedure smoothing *) 


procedure check(va 
var 


IT pty point array ); 


count , COUNTD, Councp,counth, counta integer, 
counte,countr,countv,counta : integer; 
r,C,88.,cc : amtegeae 
begin 
count :=0; 
countb :=0; 
countp :=0; 
counth :=0; 
countd :=0; 
counte :=0; 
Count r : ل‎ 
countv :=0; 
CORDE : =0; 
aa := Q; 
cc := 0; 
forr := Ito rov do 
for c: I to COMdo 
if pt(.r,c.).slope < lobound then 
aa := aa + 1 
else if (pt(.r,c.).slope = 999.0) then 
countb := CouneD 1 
else if (pt(.r,c. ).slope > upbound) then 
Gc om cocco 
else 
case VES TIC) CUT oF 
'p countp : I COURED TL: 
'd' countd := countd +1; 
hs counth := counth +1; 
Don counte := counte +1; 
E countr := countr +1; 
A countv := county Tti; 
‘a! counta := counta +1; 
end; 
writeln( level points are',aa); 


writeln('unsafe points are 2G 
writeln('for safe points: '); 
writeln('count planar -',countp); 
writeln( count hill = ,counth); 
writeln('count depression -',countd); 
writeln('count etc = counte); 
writeln('count saddle =' counts); 
writeln(' count ridge =' ,countr); 
writeln( ‘count valley = ' ,countv); 


writeln('count pass =' counta)*) 


end; 
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Main Program 


A A AAA AA A و سے‎ cape لص لحت م لومعم‎ em ت تمت سي‎ A A A AA A A A o U a O O A A AAA ce ee eee eee oo س‎ 
€ mmm ees es ee 0x 3 = 


begin 


end. 


a ei(data, data ell); 
rewrite(outslope,'data magl'); 
rewrite(outcur,'data curl'); 
rewrite(smoutel,'smdata ell'); 
rewrite(smoutslope,'smdata magl'); 
rewrite(smoutcur,'smdata curl ); 


OL 0 =S] to col do 
for r := row downto 1 do 
readln(data,points(.r,c. ).el); 
E (pones); 
121625261265 (: 


calc_slope_curvature(points); 
(* check(points) *) 


rC = ] to. col. do 
for r := row downto 1 do 
begin 


writeln(outslope,points(.r,c.).slope); 
wedpelp»ntouteur.points(.r,c.).cur); 
end; 
itor ©r:— 1 to. 10 do begin 
smoothing(points,spoints); 
calc slope curvature(spoints); 
check(spoints); 
Points := spouints; 
end; 


tor c = T to col do 

for r := row downto 1 do 

begin 
writeln smoutel spoints(.rYr,c.).el); 
writeln(smoutslope,spoints(.r,c. ).slope); 
Weitelm(smoutcur,spoints(.r,c: ). cur); 

end; 


89 


oe eee ee ee 
Sp em me mm me ec me rr rm سا‎ ee ee ee ee A ee ee ee ee ea ee عست لد‎ 


LIST OF REFERENCES 


John C. Davis and Michael J. McCullagh, Display and Analysis of Spatial Data, 
John Wilev & Sons, 1975. 


Theo Pavlidis, Algorithms for Graphics and Image Processing, Computer Science 
Press, 19387, 


Dennis Duane Poulos, Range Image Processing for Local Navigation of an Auton- 


omous Land Vehicle, Master's Thesis, Naval Postgraduate School, September 1986. 


Donald Hearn and M. Pauline Baker, Computer Graphics, Prentice-Hall, 1986. 


90 


INITIAL DISTRIBUTION LIST 
NO Opies 


Defense Technical Information Center 2 
Cameron Station 
Alexandria, VA 22304-6145 


Library. Code 0142 
Naval Postgraduate School 
Monterey, CA 93943-5002 


نم 


Chief of Naval Operations 1 
Director, Information Systems (OP-945) 

Navy Department 

Washington, D.C. 20350-2000 


Department Chairman, Code 52 ] 
Department of Computer Science 

Naval Postgraduate School 

Monterey, CA 93943-5000 


US Army Combat Developments l 
Experimentation Center (USACDEC) 

Attention: W.D. West 

Fort Ord, CA 93941 


culum Olie. Code 37 1 
Computer Technology 

Naval Postgraduate School 

Monterey, CA 93943-5000 


Euicssor veil C Rowe, Code 52Rp 
Department of Computer Science 
Naval Postgraduate School 
Monterey, CA 93943 


نم 


Professor Robert B. McGhee, Code 52MZ 2 
Department of Computer Science 

Naval Postgraduate School 

Monterey, CA 93943 


Professor Michael J. Zyda, Code 52ZK l 
Department of Computer Science 

Naval Postgraduate School 

Monterev, CA 93943 


91 


10. Captain Do Kyeong Ok 


SMC 2211 
Naval Posteraduate School 
Monterey, CA 93943 


Captain Seung Hee Yee 
Dongjak-gu Sadang-Idong 1019-10 
Seoul Korea, 150 


92 























Thesis 

TOT tee 

Gali Wipeesaleocttams for 
planar-patch terrain model- 
ing. 


VEMLU 





















o e el و و و وید‎ APN emer be deti a PELAdi d d عا‎ 5 P9 ue wie nh a oooy ب‎ 
Lia ا‎ TIAS DA RAE 2 کے‎ pps T 
a VR 9-0 APRISTA EA CNEA A P i 
2 ee | DER TT AA ANNA] AAA AS "CUPS 2d hr 2 - 
cc OSEE Ro ee algorithms for planar-patch terrai 
E ir ree d EX DO اجا‎ ANITA E AN os PTA Y PASE AA M م‎ I I ' T TT 8 K 
A ecert Pel da pie iie ا‎ a ue ا‎ P | | ti |! Pn . 
poe A Sandy som ved M a RAS ETT TETUR MNT E loda PE Pur Tu cm | | | 1 0 à 
ae ap N LE 4 fr. re (hut ef lo d orina hU Lm os om ANAIS | dad" 4A 444 [124 | | 
TT ng, PARE ILL eI Pyth go P " JEN rv TT, | | | | 
e at ET AA A ATA wae TTT PP TTT ATTE | | Il | ن‎ à j 
22 Pu AMI 1 eran 94 ARANA TA Sl | | | | t a " 
apto Porte olo dh o Ion aout ont nd t nv ey) ts eae ee De) ed s: abt est ás. | | Wil ' | - " ١ i 
eee rr eee ree ge^ we gro. HOP Tear APPS Vr gre ae ee Pe : r ELI Ii E j ; 
A A PTE A art aot s "ier : i LIBE I UE M 1 ; | 
U 8 


p, d 22 poe weyers qe egi pen | 
ILC > toed ی‎ 140 ETET ILES My ais ns me 
Tr ا‎ reg "opto Piece remi Tr^ Y e O YA T4 FJ vy ey "n 
Foe gp age. ame dapi A A Lib ate: ابسو يم‎ ere Maa 82735 6 : 
oes ^ - 8 


rip unie Tm AUTO TFNO ARTETA ETA TRA EIA A ya? 










































































































































































PO idi IS A ee eee RAT TN € PAP VP" ep 
pue: AA rrr YS my shore ERE E ARA RA PIS An pod DUDLEY KNOX i a IS any 
(sour De ee, AA A TR AA PE A OT E EN ARY ' E 
pem ARAS IAE IAS AE CACA ENS ar M ne : A da 
[nl PUTENT NU Lid E TOOL 7 e ME " swt ET ET T Mre h 1 5 we us 
AAA ATA A ال‎ PARAS TS IAE ere ser eae S| n" "TT 68 — X. 1 5 ys E z d 
pec OPH, poi pre ARPA TEA V bliss caede Ei act Peer sí 1 AA 0 a SC 5 B " cm * 5 8 s s 
prajn AAA A TTT 0 n "Pre A, ét FP - 1 E, لمعيس‎ ^ ME did | ١ 
IAS A A O * | Py br? pr pe^ Mtn A | ARS PE ae one 7 2 ni ae A s i i 
er rra pot ia PAPI ARS A TIAS AA Lees P AORTA are we A i T 5 5 38 L 
pa r ETT A ا‎ ERA. re iy ¡LARA a) e A, 1w وی ی‎ CANA AA " TI eq [NE 1 d " 1 : 
p cil heap spre ey PRE A Ee qoa ye ab hne, t FILS بام نيع ميته‎ "TI Moe wr mE PPP " ^N 5 oe 5. Us á pza 
psy Pme AA «v3 "S [E LT "YE PCR TE Vir, EXT Y A AA ٠ 1 D y" 
لسار‎ a CEEE T T AAR MD d A AA AAN AR Perr er YW PEE 5 TW 5 0 7 F - 1 5 k 
pS Capek AFRO MA METE E ore P on ms | $ ف‎ y 
A A d] 4 bo RA AT O A NA Arda ert vy hides Be rre ir vota APO AR mn beri qme i ni T ‘ E , " 5 
ا ا وشيم‎ T ا ا ا لا‎ P^ may eor ME EETRI te MARE HALE UY sA 4,8 PR! " 3 . 3 4 d 
onn iubebo aee eai A a E يلي‎ ED de Xa ARES ANA E d " P 1 hrs, o J 7 : P AE 1 $ d 
Ro e pa poop pin dep و‎ da S a an A Sa اہ ا ا ل ل‎ LA» i ety f 0 y g 
ASAS a PES EPIA + wate oM» da Stil kisa 3 4 patel 1 e 0 d A y ie d 
nane r D A A Rd agami me 1 م‎ y AT 5 3 B S Aig 8 0 5 8 
PS dede ال‎ duni ASAS A IAE” م‎ A P ^E "m " Ls 
Dr RD pee AE NY ^p 4 Je i te Pe 7 rM "et 4 " 7 A أب‎ 
A pean Ace PULIANAS AAA. لحي‎ tee cen ie , "he n - ^ isms ‘ 
lo aen pr Món PE AAA o A TEA TI NE TE 0117 . 0 , a 8 
LL EA A SÉ Unam, in DAT yee ABMS SOLED OR t سدم‎ : TP ‘s 0 š 
PA AAN T EA AIR rito a alea P » : At. z d 
mese A EPIIT" P TUE EC S LITT PPS "nnm er a y d 
P d ce aar botte . کد درجمو خي‎ aa که کم‎ aia ا‎ AA TEA AAA A aras t allt Y 3 4 
Pr a ee OT TT onl i aden bd etka sat af Marshes okt dda pa aS Pek cok ae e P tn M" n = d 
hes, ce مهار‎ ae pon AAA or TOP A Pippi) pm ARES E ا سر‎ A DP M: rr Re " A 
e a TL TT TE TI EE , ATTE IATA ISI شه | ل ال‎ WI e 9 d 1 
dd ders a) tai ntis LE Sd P PP IP T PE ACER, NS fMi T P " " P 
3 a يو‎ e< FE O AAA —— IV TM IT PIE IT O ۴ DT 7 8 zo 
PA IN O e Ped eet teeta es ur TL Pee لس‎ "M" 1 > P i #7 
p oor or rm ne OTT ee ss a EOI EP eee) oe ey Tor E eee E pe pa J 1 " 
A qui cres A AS EE OR NO TR ES TEN WP 9 " x " " 
ines atii iri sive poy oe yore. DAA IO ET A apice PB oh rer ان لديا‎ G 4 " FU هم‎ * 0 ^ 
Hee orar ibit raid aL rati in Im PETRUM CC m T TARTA A Pm "E ١ T " 2 , 
رورسم‎ AA A ا‎ o A ep TS aaa UT aO Lu Dui iss A A] E , A P y $ E 
ELA a OR Pags AI AO و پک‎ CIA do AAA ti A get AA AA G LXI "TI" ^ ' : 
pr insti A e Er TEE RP rer epee rer T aA ATREA o ME Mir" "un > : P s 
A AR ا‎ eRe AAA ART AAA AN PARIS AA Mala tof as 9 5 teg 2 E n 01 T 
Ds pr ares APARATO ET a deed IAE A E T 5 "EI ' » 4 44 8 
PA AA ا‎ d a A eh O TTT eT OTT IA ATA 3 omen 7 0 m TIT "IT 1 4 g x » 8 
ومع ضرم رصن نوصح‎ AAA A A م هس 4م فش يوبن‎ 147 PA OEA PRO AA ot! * ah ‘tt? Pepi er ee y P f 8 Pa a 
RAI Ed A A A ل 0 ا‎ 0 e a ETT E er T 1 , "um 1 
paca bat pt Sage E e ب عر‎ PC zm A OPT: 5 : 1 i 8 A . 
A oso It at nt P ao vl d unt harten Mae dh me oid eB hun s oim ¡ERIN e DP ^ pe eee mot CN» ps er SES h " PARLA 
Peon e dio نا حي سارو او‎ E AT ا‎ > Fer ER Heo - - $ j A O 0 
AAA A A RO "vi VADE PRU E PP of imer MSAN ^o^ "EP UE T AR i d = 
E A AA dise aid bt dt dn RUE a E ITE utri ‘ P PERTE 5 , Ne Fae i 
PAS da P mad onde Lu M MULT id Hoy ERITREA TI AN M o Aw! 4 "n "a 1 8 5 
poser AAA AT EE E Pm D e4 d AID ll 0 write A fap 1“ # E £ a d 
adf o fà uoo dk oh f um arme e go og oco en nae pao mns Y tT Er UT 0 ir E ELA "MIX i v2 r 0 a ates E ALL مايه‎ aa 5 A it ke ^ 
d Td ode سخ صر‎ repite era eather oan oa AS ie M ou y Uoc mE ilf cue Noo TIPP PED ^ 4 ^ E ; 1 : i 
pi خا‎ el di EE, den A م‎ A posi SE PRSE ^ i "ann ‘ys 1 f ovn AA eur T ETE 1 d a 
red PPP d PPP ITD OTT TT TCC OID esi dL ERSTES amass PTT dep rt o bury, "ERT BEII AD 18 4 pda pr IAS 4 n 5 4 f 5 ^ 
ل‎ omae al vier lin pm TP Men tae! IA ce ee y M 9^ AMT LE TP PT Lo 4 " 254 31 TO Pd P 8 
Pr IA IS TIE TA NS UM TS TI ÍA PANA y m "SPP و‎ on di F S. á 
mgt wo pp ga pice cape dns byes PO FPG Ls SS POR Sa Nea Te TF as Amin EA AAA het oh k Jas TE M Cre TAR a hp , " y DS 
d 72 0 الل‎ A A PATIO O Pa m E - "mov prr [ ^ PI ITE ee) 2 r "TI NAM 0 "P 1 es. E Py 
ds AO ya HP prr ee ETE re TIAS ET A Oe ae 5 "e TU 3 "WP ee ^ d i 
0-0 p te gr wee DESI AAA 0 AE PEE Or AS e Y Pe poc ا‎ gg. , &4 "mew re T. 5 3 d 
Lu vel ei Bie eT al ole hacen tetera JA P TP يتا وه‎ PON ا‎ ALA e. 5 0 1 1 2 
Ei A A pita penne ت کار‎ eds ls e او‎ ARI EET TP dec ar as lar > "a " de 9 > 2 y ee 
oda pá ba o rr or * عاب‎ ٠١ صر في ىر‎ pe an PM On اعم ب‎ APD A VERD PER "يح‎ & Mf oer Eoo A B ERU NOE 3 j 4 3 ne 9 
¡A e mer — m peere aT AME TC TP our 1 T 2 JUN $1*9 sg 5 ae "TET » "PIT es ze f : d 
لون‎ 5 gf LÀ: Fut amid ES E pe Pine, a) ihe ha $ Mal od Lai nin, LITT M " rin At A "7 oe Ne T Sd > 
sr” PLA ا‎ eet sic: of JA ^ o ARA Qu T Ei ا‎ LE m DIETE Pa L] y s 97 P " 5 , e 0 9 A j 5 i - 
EDGE A ل‎ FE id dz KE La ae PL QV P PT le ms O ore PV E fd " TP yr ed die od ٠ ` 3 
m: PET TEET al da! y alr موم ار‎ AAA TE ARA er A ey AS "af yt 3 x d d d a 
D) o IP Ped b ind ad Seti c e o m dO *» I9h5 AP i م‎ ETIT E e 5 D a 7 Se Hy 4 5 - 2 E , , 
mesa r ry Li grin. Toe ow yo. اود کې و‎ m ¡PAT m pr e Py "m "TE " 3 An x 
DS ae + "pb > as T MS 3 أيه © 4 ريت بهد‎ L] PRSETER i 4 dt 0 AA ا‎ "mo i 0 
9207 ger a a HR A "el gf. we ' اص‎ "i pete ا‎ LP" ay be E "POMPA PAE A ^ni K t . 0 0 0 P 4 0 
AA AAA ا‎ De A LLLI :م»‎ LIC zeda IP m 8 A TE - رك 2 به‎ E ا‎ 4 Purr 5 0 y 
M o Tao ll o a جه‎ Fe =~ d Ra dil dod A nm PT e" nig A s 0 à o 8 
ا ا‎ Spal tk ebm AAA مو‎ gir ZIM PIU PP E ILI CY PORET SN fact TO 8 0 1 5 K j i 
Del ro hr ee RR "nee v rmt. ^ hb 4 "TP Ma Pers PAE 2° LP "m rir dy b GPS v PT 1 " a 8 er. 1 "eum 
dE E "du ww ¿e dy a 2 CAY j RS ae LS LAN 45 d d O 5 t d ` 
de Joao Ae arr MEX erated TA we a RA PITE "e 4t 4 (hoe Ar 3 t E MEA E 
"a عد الها لوص آل‎ br a > h A ete £o a Fur PY ot a ee oa ا‎ "m "TTE PT "I. " A $ 5 
ah Lo A E Dj e I EE i lk E A y 5 2 1 t qtio j 7 ab D $ 0 
ا ا لا ل يدت‎ CC Y d 
Eon dr a LEIA 4 = d^ s E se 8 y 5 
ت حر ا‎ o de ds ad E - 0 0 d à & si e i 3 A 
d ا بيع‎ TE FS Wi 1 Ea un 
ہے‎ EIA 3 A N 1 
E or as tery ana 7 8 s 
EL OEE es meee ' ^ 0 i ' D A " 
x oa she Lacs ow عه ل‎ A s A E Pl Á 4 0 ow : 
= ms 3 oa, d 
- —— T3: Tr rectam "d pr i a E 
ore PP a nr A n. AL ١ y id 
e ٠. في‎ * 7 tne LI P] ل‎ 
, LI LI P n , 01 0 
L] L] * LI 
« g 8 5 E 0 , 9 
1 L] A L] LI 0 T a 
= E A O هه‎ 
n ^ EE i 4 i 
Ito ae) a 
d 8 "a" 3 0 8 "E 
" d a "T 7 
, 31 5 8 1 LI " 
3 as 3 °" م‎ 8 
Gd L] 
4 لي‎ 
LE : " 1 0 
L 0 0 à 3 à, 2 ے‎ 0 
E ey P 8 0 " J 
4 E 0 7 
D 0 d P 5 ô 
o $ A 
o. 5 : 
m 7, 
Eo, 8 8 
, ql , 5 r 
8 " " 
3 ¢ g ' "En 
8 
7 > Di ^ 5 8 
» 5 
Li , 0 a ° 
e , 
Li LI j E r 
1 0 D" ri Ld 
$ tos 
"I PV n P 
ee 5 " 7 a 3 
pe L] 
2 P . 0 
L] MI ل‎ RE 1 d 
d 0 D e 
ا ا‎ E i 8 P 
Pe 3 3 0 , " 
, € Y ۴ 
1 4. L M 
"FLUR no E] 0 $ 
r1 E A 5 
S X j 1 
5 e 4 
" " 0 ل‎ 0 
LI ام‎ 
3 8 
A Ave " ^ , 7 
L M . 0 1 5 
P 8 A F 
7 
' 
LP 5 
tal, LI d 
i Ti 0 wu 0 
e) "n E B 
Cur dd y z y ۰ 
' d + ns 
A ^ 
, " 
^ 0 
be 4 d 
q. d. Toi? 
dela e EP M P 5 ied: acra 4 M 
A e 54 s. aeri 8 Ly y as ١ 
| ونا‎ a e y ra q Sa (o tel PA rud ١ 
a dy e T pe, had t E D P 4 
eee ا ا‎ AS te يه كه‎ j 2 
q co tl Mr e o مدن‎ Uer ho dl Lae - چ‎ f F : 
a ie SISSE fe tw ey TT ew + 
acp IT عه‎ i ve^ “ 58 
"f OPNS EE A ML LOU Pr et ١ , . o i 
$v DANE ا ا‎ wut 
ION ro try kia xa D 5 3 0 £ 
ا‎ od “al Ei a: < E ira d T è E 
immi Po twi EN AS vant y ا‎ + 4 1 3 : 3 i , 
IDEE E t SA PO iodo teh Oh aaa, " Mi PAL! N دان‎ 
3 rte n my ate gyra YAA ل ا‎ 0 den M D 
TE eR e poes 1ن مرا‎ id lege mie Hwy e : 
c mae t am qn e nanan da d oM. das a Lada dos M "oin 1 
ea LS ree sa be ahh ae ay get t e: ey 2 Lan i ; 5 
eta cn bab A م بج‎ RES A pc LIE ا‎ F A ~ d 
ا‎ E ا سن‎ j aiaia qe Vr : e ET uc 
ORUGA 1 and £e. a 5 
a on ata cote Del ej Steak cr ees 1 " t M ; 
4 3 Ld 5 
a AT ا‎ anh Tae ahd tas Ar n H A 2 
rm o Se Téc Di يه‎ a re d wv WA LT n Pt F ١ 3 9 
lite te nel ao M. al! mpeg Lê kê vryt yur A a ؛‎ * Y* m : L 
EET! er alg ee zi a pepo tele aie a Form vars 4 + i ١ , y 
A ds تروب مج وس‎ da Ba ھی وا‎ es M 35 0 4 ٠ n g » 
7 ct ap ee Op tll pue aca ie Ma TO RA وام‎ ie PAS "ea Le iad be wey 3 5 , 5 EA 
vie Sy og "a n^ ander st hip gares e i c, Pri baci LL اله بدي ريه‎ ry : ‘ 4 y 
ITALIA Y? A Ano) PRA AR A Ao ur rt $104 3*8, V Mut (e .t E wo! e ^ = = Ht! : 
Mu sk bid im Mei Là i Mrd de ye n , 1 ۹ "ECT . "Uc we d : 
1 lel ed ard. ata ire ^ Aie) Aja el p IA dd dan Aa ١1ج‎ dh A | 0 E 
La] 4 - SA den AN Li alil. ي‎ M *a vfi vm Fi Pay 2 a ea aiia, i 3 eu 2 j 
epu empre de rtm Lone ron yt atados a wf ety IR ١ © pá M N: J ر کیک‎ ws $ E Iun ms “y 4 
دسل‎ Lais mea ac oye vein ett SIM AL T ay us PL Ath etd RU v5 DA T Ei A A ا‎ 
KS Cot poe i ue ie AA id "A nat atom fa aes ids ER Po [ , 8 " me d E 
orate be ee 3 LARA" ha sa Uadec de %, uev tervi I ; b 1 i M 
J "b 7 Tus p LR 04 C La a 1 «i d 
arty by + A ل پا‎ n rd x d nf 0 AOS k^ Ardo 3 i i E y ; i 
ae LEE SUR ari oe Pcr Tc NAA It +4 lg ra 8 n i 
ae aed ^ ues ^ ee poe e e e iat ré es Aare eee ys ا نهنا‎ © ١ ‘ 
- b A vv mE. 
DA e are as URL) E EA Yu e EIC HERREN d ١١ 1 ii 
E LA A A Coh EET AIB AU M OA ER at d 
rv a fe ty kat rey ساسا‎ wal the ta aa bred sadi Ye ¥ ets UTA fe bo ١ TO ia ^t ل‎ y = d 
0 AAA Abd A E lere Ares eii OS CIL IU A EL NS We v e " 5 
oe tle le eb E A E ee) Aga LAN ie te 1 TA SL A L d 59 8 
مضت‎ Mi Re A vana I e rrr m adeo d toot do n Reeth Sed a a ETEA EE s o | -n y 1 
ei * ee Pi هرج‎ ña Y IA TOT a NU " "m Be م‎ Y " 
A PU و‎ A oe eer lone es yr ere ye tei heen red di Bed v 4 E ۶ E a had 8 1 4 D s 
in iaa mda eR hii Sae dà ose dn N ia AA NE ene jy ey ts who e > 1 au S , ee b 8 
ee od ca le hp bh ahh kh ences Ay ا‎ wl Sots DAE E RA 7 A RACE: RAS ١ : i , 
Rd ا‎ be alec y Lg tan wus Dh M" MOL YA Ct يان‎ RA - y vs A : 
aw یی و ج چچ‎ uii ١ TETT AES 3 e rs th D ete aD qe “eb ata. Wt . "T Y a 
Phat Sd كيك وسو بودي ب‎ co sev ts 3 SS de A ae e vt A 6 NI E 
E ACTU» RAI LE لفن‎ y Pro, NY Pu dA A da Q4 e Uer 5 5 35 3 1 i e 
mL PRA t cido pla A (er I Y-€ ١ 7 ‘ 
pr LA Lied ا‎ A N Bd EOS EA A Li AA uA D ^ t 1 
CER pode LL arde d مضه‎ y ip] ns re dy" e be d کا‎ 2 x | 
"m AAA tr kes, E toads eder t M oL 0 g 
a عه‎ bis, A p > A N dridi A n vro vr 0 É 
RE yd ra ej Da al SALEEN s > pM " ! Gv" 
al nom] AAA العم و‎ ۹ A Y o o A = 4 i 1 ; 5 ^ 
om m due n a oia La Aid do T dicta ña et teh A alia LI 4° "M AA q. NE err Y ١ 
DER ا‎ La PONSA dl id M D i aD g M 0 > AA Y kl 
te کے‎ we Ah 4 n Eo pa CETINA Y 0 P AER . 1 هه‎ 7 n^ L 1 1 
RI tiden de Ae irae f» 9 y 7 n" a vr Ye | . i 
T! ERAT ce eA R vr? A 2 5 
SED AG pada ga + ee hid oka hte a S ¬ ` 3 a! 1 y 
d Loteri Amei Wa ahr dl wy a A! ” m Se * ve 3 K Š 4 E 
O WS: IE Pies 4 عر‎ Li tes ie Ea s 
a E 2 x ile 
> Ae pe TY we E P : M À Mir: . A i 
"$ ^w» crest V 4 » 9 y 3-9 9 - 
M p vieta m or’ 2) ا‎ V M 
Dele MEE Ex 4 - ?, *y LA d y 
i ent E dc el “٦۹ 06 Eb» ١ 1 
لط وما اه‎ A AA a NA y rt L x 
mine d t Y é 7 


anyen‏ ناه 04 بىر DS‏ شد م 


E M ~~, A e p " F 5 
PENE ERIN A A "y ev o i 5 


